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lf INTRODUCTION 
A. BACKGROUND 


The passive geolocation of electromagnetic (EM) emitters plays an increasingly 
important role in all three aspects of Information Warfare (IW). In order to protect 
friendly communications from hostile jamming and interference, information warriors 
must first locate the source of that jamming and interference. The geolocation of an 
enemy radar exploits the enemy by gathering information ascertained from his own 
information gathering system. Finally, in order to attack an enemy communications node, 
its coordinates must be passed to the strike aircraft or the cruise missile. In addition to its 
importance in the military arena of IW, passive geolocation of EM emitters finds use in 
law enforcement surveillance, search and rescue operations, navigation and the 
enforcement of communications regulations. 

Three aspects of EM waves lend themselves to exploitation for geolocation 
purposes. The path of an EM wave can be closely approximately from transmitter to 
receiver given the frequency of the transmission and some basic characteristics about the 
atmosphere through which it traveled. In addition, when EM waves arrive at two moving, 
spatially separated receivers, the receivers measure different Doppler shifts in the 
frequency of the transmission. Finally, the time that an EM wave arrived at two spatially 
separate receivers yields valuable information. These three characteristics naturally 
enough, are the basis for the three basic methods of geolocation, Angle of Arrival (AOA), 
Frequency Difference of Arrival (FDOA) and Time Difference of Arrival (TDOA). 

AOA geolocation involves the use of highly directional antenna arrays. By 
measuring the phase difference in the EM wave at each antenna element the receiver can 
calculate the direction from which the wave emanated. Drawing a line from the precise 
position of the receiver in this direction yields a vector containing the position of the 
transmitter called a line of bearing (LOB). Moving the receiver and taking another 
measurement or using an additional receiver located elsewhere produces another LOB. 
At the intersection the these LOB(s) is the transmitter. AOA is the oldest of the three 


methods and has been known by many names throughout the years of its use. Beginning 


just prior to World War II, the name most commonly associated with AOA was 
triangulation. With small measurement errors, the three vectors had a finite width when 
drawn and thus intersected to form a small triangle within which the emitter was located. 


This point is illustrated in Figure 1-1 below. 
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Figure 1-1 - Angle of Arrival geolocation 


For moving collectors, the Doppler effect ensures that the receiver will measure a 
different frequency of arrival than it would as a stationary collector. Given two collectors 
that measure two different frequencies of arrival (FOA), the difference of these FOAs 
(FDOA) can be used to determine the position of the emitter. FDOA produces a locus of 
points called an isodop along which the emitter lies. This isodop represents all the 
locations from which the EM wave could emanate and produce the difference in Doppler 
shifts measured between the two receivers. Figure 1-2 shows a two dimensional 
depiction of a family of FDOA contours in the plane of the paper. Each contour 
represents a specific FDOA between the observers. The emitter could lie anywhere along 


the appropriate contour, and thus multiple measurements are required to determine its 


location precisely. While these contours can be approximated as two dimensional when 


the emitter is on the surface of the earth, for an emitter in three dimensions they would be 


observers 


Observer 
Motion 





Figure 1-2 - FDOA contours 


surfaces. A receiver wishing to use Doppler for geolocation purposes must possess the 
ability to measure frequency more precisely than the smallest Doppler shift expected. 
Otherwise, the small Doppler shifts go undetected as they are below the noise in the 
frequency measurement. This has been the limiting factor in the usage of FDOA for 
precise geolocation. Measurement of frequency precisely is much more difficult and 
costly than the measurement of time precisely. 

The precise measurement of time has been revolutionized by the high speed 
digital technology of the Global Positioning System (GPS) and synchronization circuits 
for atomic time standards. Using the difference in arrival times of an EM wave at two 


spatially separated receivers, a locus of points called an isochron includes all possible 


locations from which the EM wave could emanate and arrive at the two receivers at the 


two different times measured. In two dimensions, like the situation that approximately 


s(t) + n,(t) 


Observer #1 


Estimate 


TDOA Estimator of 


t-D)+n,(t 
Observer #2 s( ) + alt) 





Figure ]-3 - Generic TDOA estimator 


exists when an emitter is constrained to the surface of the earth, the isochron can 
be represented by a hyperbola. In three dimensions, the hyperbola becomes a hyperboloid 
of revolution. A simple block diagram for a generic TDOA estimator appears in Figure 
1-3. 

The traditional methods used for TDOA determination include the general cross 
correlation (GCC) (in the absence of Doppler shift) and the cross ambiguity function 
(CAF) (both TDOA and FDOA simultaneously). These two methods perform well in 
benign signal environments. That is, with a signal of interest (SOI) with a high signal to 


noise ratio (SNR), in the absence of jammers and other signals not of interest (SNOJ) in 


the frequency band around the signal of interest (SOI), the GCC and CAF methods are 
superb estimators of TDOA. However, in a peacetime EM spectrum increasingly packed 
with signals, and a battlespace spectrum potentially flooded with SNOI(s), hostile 
jamming and weak SOI(s) the GCC and CAF methods suffer from severely decreased 
signal to noise and interference ratio (SNIR). In light of these weaknesses in hostile 
environments, a relatively new concept in signal processing called cyclostationarity has 
been applied to the TDOA problem. 

Cyclostationary processes are characterized by time-periodic statistics. A wide- 
sense cyclostationary signal specifically has a mean and an autocorrelation that vary 
periodically in time. These periodicities are usually hidden in traditional second order 
stationary processing techniques such as the autocorrelation function or its frequency 
domain counterpart, the spectral density function. However, when second order 
cyclostationary techniques such as the cyclic autocorrelation function and the cyclic 
spectral density function are used, wide sense cyclostationary signals readily display these 
hidden periodicities. These periodicities now revealed may be utilized for the detection, 
classification and parameter estimation of the signal. Among those parameters that may 
be estimated is the TDOA of a signal impinging upon two spatially separated antennas. 

The robustness of the cyclostationary technique for TDOA lies aan the fact that 
every signal has a unique set of periodicities that depend on such parameters as 
modulation scheme, bit rates and spreading codes. Two signals that may appear as 
spectrally overlapping in the traditional stationary spectral density function can be 
separated with great accuracy by the cyclic autocorrelation function. In highly corrupt 
environments, cyclostationary processing techniques provide the capability to select 
specific signals for geolocation, nearly independent of the presence of jamming, inference 
and SNOI(s) that may spectrally overlap the SOI. 

The United States Navy has historically relied upon AOA techniques for direction 
finding. Beginning in the 1950’s, the Navy began constructing the High Frequency 
Direction Finding (HFDF) network of circularly disposed antenna arrays (CDAA). These 
CDAA(s) consist of a circular array of elements and reflectors that serve to detect the 


direction of incoming energy. Given the mass of recent base closures, their numbers have 


dwindled significantly. In addition, they are not an organic asset for battle group 
commanders considering the current tactical communications environment that relies less 
on long-haul HF and more on Very High Frequency (VHF) and Ultra High 
Frequency(UHF) point-to-point links and Super High Frequency (SHF) satellite 
communications links. 

In addition to the HFDF network, the Navy added the OUTBOARD and 
OUTBOARD II systems to some of its Spruance-class destroyers. These too are 
primarily HFDF assets and have been followed by the current construction of COMBAT 
DF equipped Essex-class amphibious ships and Arleigh Burke-class guided missile 
destroyers. These three shipboard systems have proven extremely useful during the past 
years but as they have become the object of intense scrutiny with shrinking ship-building 
budgets and increasing decommissionings. Consequently, the Naval Security Group 
Support Activity (NSGSA), in June 1993, contracted with the Applied Research Lab at 
the University of Texas at Austin (ARL:UT) to develop an affordable, low-risk 
TDOA/FDOA geolocation system using commercial off the shelf (COTS) and 
government off the shelf (GOTS) technologies. 

The Carry-on Multi-platform Global Positioning System (GPS) Assisted TDOA 
System was tested for the first time in August 1995 after fifteen months of research and 
development. In further testing off the coast of San Diego, California, the prototype 
system geolocated HF, VHF and UHEF emitters with a mean square error of approximately 
100 meters, using a CAF based TDOA determination algorithm. 

As noted above, the CAF has performed exceptionally well in the relatively 
favorable conditions during prototype testing but in theory, is susceptible to low SNR 
conditions and co-channel interference. Cyclostationary processing techniques for TDOA 
determination offer potential improvement to the performance of the ARL:UT system in 
highly corrupt environments. The evaluation of a cyclostationary TDOA algorithm in 
conjunction with the closed form geolocation algorithm used in the ARL:UT prototype is 
the primary objective of the work presented here. 

As a secondary objective, the cyclostationary TDOA determination algorithm and 


the closed form geolocation algorithm will also serve as the first two blocks for testing in 


a developing geolocation software workbench at Naval Postgraduate School. The intent 
is to model the system in the Simulink® environment developed by The Mathworks Inc. 
and allow different algorithms to be exchanged by merely interchanging the appropriate 
Simulink® blocks. This will allow for the testing of all aspects of the geolocation 
problem from data conversion, filtering, AOA, FDOA and TDOA determination and 


various geolocation solutions. 


B. ORGANIZATION 


The traditional TDOA determination algorithms, the GCC and the CAF along 
with their disadvantages in corrupt environments are present in Chapter II. Following the 
illustration of the vulnerabilities of these two algorithms, the theory of cyclostationarity 
.-and TDOA determination by cyclostationary techniques are presented in Chapter III. 
Chapter IV depicts the closed form solution to the TDOA geolocation problem developed 
by Dr. Petre Rusu at ARL:UT for the prototype geolocation system which is described in 
Chapter V. 

Chapters VI, VII, VIII and IX constitute the bulk of the investigation, 
summarizing the new algorithms developed in MATLAB®, the test plan and the results 
of those tests including the plans for Simulink® block development. Finally, Chapter X 
draws some conclusions, offers explanations for the anomalies encountered and suggests 


areas for further research. 





HI: TRADITIONAL TDOA DETERMINATION 


A. GENERALIZED MODEL 


For the development of the two traditional methods of TDOA determination, 
assume a Stationary signal from a remote emitter impinges upon two spatially separated 
antenna elements as illustrated in Figure 2-1. In the presence of additive white Gaussian 


noise (AWGN) noise, the two signals at receivers | and 2 can be modeled as 


x(t) = s(t) +n,(t) (2-1) 


and 


ae) vi (2) (2-2) 


where n,(t) and n2(t) are assumed to contain only AWGN with no strong SNO\K(s) in the 
frequency band of interest for the SOI, s(t). A represents the complex valued relative 
magnitude and phase mismatch between the receivers; D is the TDOA of the SOI 
between the signals. 


It follows that the autocorrelation and cross-correlation functions are given by 


R(t) =R,(t) +R, (7) (2-3) 
R,(t) =|A| R,(7) + R,, (7) (2-4) 
R,,(t) =A-R(t-D)+ hone (Ct) (2-5) 


and the respective spectral density functions are 


Ome, Cy) (2-6) 


S,(f) =|AlS,f)+S,, (A) (2-7) 


S, fI=A-S,(fy-e = +S, GF) (2-8) 


These relationships contain the parameter of interest, specifically, D, the TDOA. 
The next two sections show two different methods used to extract D from these 
equations. The first, the General Cross Correlation (GCC) function, can be used in the 
absence of relative motion between the two receivers and the emitter. The second, the 
Complex Ambiguity Function (CAF), can be used to estimate FDOA and TDOA 


simultaneously as is necessary when relative motion between the receivers and emitter 


eX1Sts. 


Observer #1 
TDOA, D,, = t, - t, can be 
used to find ®,the angle of 
arrival of the EM wave wit 
respect to the baseline rj, 


Emitter 
Z 


Path difference Ir,, - r..! causes time difference of arrival Observer #2 
for EM wave measured between observers | and 2 





Figure 2-1 - Signal model scenario 


Of key importance to the development illustrated here 1s the statistical 
independence of s(t), 2;(f) and n2(t) and the absence of in band interference. While it may 
be argued that the signals and noise measured at the two receivers can be correlated to 


some extent, that scenario greatly complicates this case and is beyond the scope of this 


10 


demonstration. The derivations that follow are drawn from [1], [2], which treat the GCC 
more specifically and [3] and [4] which address the CAF more so than the GCC. All 


assume the three components to be uncorrelated. 


B. GENERALIZED CROSS CORRELATION 


Applying (2-5) in the absence of SNOI in the SOI band of interest and statistically 
independent s(t), 7;(t) and n2(t) yields 


Po) eer DK, (7) (2-9) 


This function will peak at T= D, the TDOA between receivers | and 2. Because n,(t) and 
n2(t) contain merely AWGN, their cross-correlation term in (2-9) above reduces the SNR 
of the measurement but does not add interference in the form of SNOI(s). Further 


analysis using (2-6) - (2-8) in a similar manner reveals 


OO an) (2-10) 
Ser) Al S| tS, Ch) (2-11) 
SO DewAG (ce Mat SCL) (2-12) 


To extract D from (2-12), first note 


B B 
S.(f) = a Fo A =I Sto + (2-13) 


='() otherwise 


which assumes the SOI effectively occupies a finite bandwidth B around the carrier or 


intermediate frequency fy. Taking the ratio of the spectra and using (2-10) with (2-13), 


1] 


See A-e 9? f-SSf <f,+> 


(2-14) 
SF) 0 otherwise 
and taking the inverse Fourier transform of this ratio gives 
for 
S, | 
(Ce) = | Su JF) et Ph df 
_B Sf) 
fo, 
(2-15) 
~~ Asin[fn B(t-D 
= Bada hacia cos[2 nf, (t= D)| 
m(t— D)/2 
which peaks at T= D. This in turn may be rewritten as 
lene 
nt) = JW(f)S(f) af (2-16) 
uP 


The weighting function, W(f), is defined in this instance as 1/S,(f) in (2-16) which 
is the best choice given no a priori information about the SOI [5]. In addition, it 
distinguishes this case as the generalized cross correlation method. Assigning a 
weighting function, W(f) = 1 reduces the TDOA determination to a simple cross- 
correlation as in (2-5). Given prior knowledge of the noise and interference 
characteristics, other choices for W(f) include the Roth impulse, SCOT and PHAT among 
many others [3] which are designed to reduce specific noise and interference problems. 

It is clear that, in the absence of significant noise, co-channel interference and 
relative motion between receivers and emitter, the GCC produces the desired estimate of 
the TDOA of a signal from the ratio of the estimates of the spectral density function of 
the signal at one receiver and the cross spectral density function of the two received 
signals. Figure 2-2 above illustrates the generalized cross correlation process. The two 


filters, H,(f) and H2(f) are specifically designed to remove out of band interference. 


| 


While in practice the GCC is often done in the frequency domain, Figure 2-2 depicts the 
time domain correlation of the input signals. The processes are theoretically equivalent. 

Clearly, the presence of an interference signal in the bandwidth B defining the 
spectral densities would corrupt the estimate. In addition, as is shown in the next section, 
Doppler shifts at each receiver also must be accounted for in order to accurately 


determine TDOA. 


Maximize over 
eG 





Estimate of TDOA 


Figure 2-2 - GCC block diagram 


C. COMPLEX AMBIGUITY FUNCTION 


The complex ambiguity function (CAF) can be interpreted as an extension of the 
GCC for moving transmitters and/or receivers. Stein in [4], showed that in order to 
properly determine the TDOA between two receivers in the presence of a Doppler 
difference at each receiver, the spectrum of one of the signals first must be translated in 


frequency by an amount f equal to the Doppler difference (FDOA) measured between the 
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observers. In order to show this, a Doppler shift is introduced into the generalized model 


from above. Equations (2-1) and (2-2) can now be written as 
x(t) = s(t) +n, (t) (2-17) 
y(t) = As(t - D)e"?"" +n,(t) (2-18) 


where fg is the Doppler difference as measured between observer | and 2. 
Now, in place of calculating TDOA with the two dimensional GCC from (2-16), 
which in the presence of significant Doppler could peak at a value that does not truly 


correspond to the TDOA, it is necessary to calculate the three dimensional CAF given by 


T 


A(D, f,) = fx y (t—-D)e?™ dt (2-19) 


0 


The simultaneous determination of the TDOA D and the Doppler shift fa causes 


IA(D, f,)| to peak. At fa = 0, the CAF reduces to a GCC problem as outlined above. For 


fa# 0, the CAF may be thought of as a GCC performed after frequency shifting the 
spectrum of y(t) up or down as necessary by an amount equal to fy. A block diagram of 
the CAF operation appears in Figure 2-3. Note the similarity between it and the GCC 
diagram of Figure 2-1. 

The three-dimensional width of the correlation lobe 1s directly proportional to the 
accuracy of the estimates of TDOA and FDOA. Stein further points out in [4] that the 
variance of the estimate for each parameter can be related to the noise bandwidth, the 


signal bandwidth, the integration time and the effective input signal to noise ratio as 


follows 


Sioa => ——= (2-20) 
5 = 


Maximize over 
T and f, jointly 


Estimate of TDOA 





Figure 2-3 - CAF block diagram 


1 1 

Oo a (2-21) 

FDOA is J BTy 

where 
B = noise bandwidth at the input of both receivers 
2 Ys 

| £°W.(Aaf 

8 (2-22) 


— 
[w.paf 
where W,(f) is the spectral density of the signal as shaped by the receiver. 


7. = rms integration time 
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tatty | 0-23) 
12 i Sieaaie 


I ] 
where — and — are the SNR(s) for each receiver respectively. 
Yi We 


D. PERFORMANCE 


Given (2-20) and (2-21), it is clear that both the GCC and the CAF can be 
rendered ineffective by moderate in-band noise. Figure 2-4 below is an illustration of the 
theoretical standard deviation of a TDOA measurement made at various SNR(s) and 
integration times. It considers only the errors introduced by AWGN according to the 
theory presented above and includes no error from any other source such as 
instrumentation or machine precision for example. Because theory dictates these values 
as the minimum errors, they may be considered the lower bound for this scenario. 

The model for this demonstration assumes a rectangular signal spectrum with 
rectangular, full bandwidth Gaussian noise with no Doppler between observers. It 
assumes no magnitude or phase mismatch between receivers and the same SNR at both 
receivers. It approaches the ideal situation given the perfect match between signal 
bandwidth, noise bandwidth and receiver measurements. In reality, the signal bandwidth 
would be more Gaussian or raised cosine pulse shaped and the receivers’ measurements 
would not be matched in phase or magnitude. These differences would corrupt the 
TDOA estimation further, requiring more SNR and integration time to achieve the lower 
bounds determined by this idealized model. 

Clearly, below approximately 10 dB SNR and 400 ms the idealized model 
exceeds a TDOA standard deviation of 100 m. Given three TDOA measurements at this 
level and a simple Euclidean norm combination of these errors, a geolocation derived in 
the absence of any geometric dilution (see Chapter V) would have a standard deviation on 
the order of 173 m. This exceeds the 100 m goal of the ARL:UT project without 


considering the additional errors introduced in an actual system (again see Chapter V). 
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Figure 2-4 - Theoretical standard deviation of TDOA measurement given various 
SNR(s) and integration times 

Poor SNR(s) and SNIR(s) are obviously problematic for the GCC and CAF given 
the need for geolocation accuracy on the order of 100 m. While increased integration 
times can solve some low SNR problems, this solution has an upper bound. The 
integration time must not exceed the coherence of the SOI. That is, the statistics of the 
SOI must remain relatively stationary during the integration time for the GCC and CAF to 
operate properly. A significant increase in integration time to combat poor SNR may 
exceed the coherence time of the SOI and introduce yet more error. Modeling the SOI as 
cyclostationary vice stationary and employing the appropriate processing techniques can 
in many cases overcome most of these problems. 

Cyclostationary techniques exploit periodicities introduced to man-made signals 
in a number of ways. These periodicities can unique to specific signals and thus can be 
used to distinguish one signal from another in the same bandwidth as well as significantly 


reduce the level of post-processing noise. Thus, with cyclostationary signal processing, it 


Ly 


is possible to tolerate significantly lower SNR(s) and still obtain excellent TDOA 
measurements. Chapter III introduces the theory of cyclostationarity in communications 


signals and develops a method for TDOA determination. 
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Il. CYCLOSTATIONARY TDOA DETERMINATION 
A. DEFINITIONS 


The concepts of cyclostationarity have been examined in theory for over two 
decades. Beginning in the late 1960’s, Dr. Franks of the University of Massachusetts at 
Amhurst and Dr. Gardner of the University of California at Davis (UCD) began extensive . 
research in the area of cyclostationary signal processing. Dr. Gardner has since produced 
multiple publications in the field. His paper on the general theory of cyclostationary 
signal processing, published in the April 1991 edition of the IEEE Signal Processing 
Magazine [6], stands as the key reference for most aspects of the theory. As is always the 
case when one individual has contributed so much to an important engineering 
development, much of the theory and examples presented here follow Dr. Gardner’s work 
closely. His solo efforts and collaboration with others appear in references [6], [7], [8] 
and [9] and is redundant in many instances. Where appropriate, the specific source of 
unique information is referenced below. 

As previously noted, most modern signal processing techniques associated with 
communications applications treat SOI(s) as stationary random processes. However, 
because most manmade signals are generated through some repetitive, periodic process 
such as the amplitude, frequency or phase modulation of a sinusoidal carrier, the 
encoding of data or the encryption of a message, their statistics inevitably vary 
periodically with time. While in many instances, receivers may successfully ignore the 
underlying periodicity of a manmade signal, often the detection of the signal and the 
estimation of its parameters is more successfully accomplished by modeling the signal as 
cyclostationary vice stationary. 

Simply stated, a process with statistics that vary periodically with time is termed 
cyclostationary. Figure 3-1 depicts a block diagram of the procedure that leads to a 
cyclostationary signal for most communications processes. A stationary random message 
such as digital data or analog voice is modulated, clocked or framed by a periodic 
communications process through some nonlinear system and a cyclostationary signal 


results. 


Ae 


Stationary 
Random Message 
- digital data 
- thermal noise 
- analog voice 


a 


Nonlinear Cyclostationary 


se System Signal 


Source of Periodicity 
- RF oscillator 
- bit-stream clock 
- repeating frame structure 
- rotating machinery waveforms 
(motors, engines, turbines, propellers...) 





Figure 3-1 - Origins of cyclostationarity, adapted from [10] 

From a strictly mathematical point of view, a cyclostationary signal of order n is 
one that will have additive sine-wave components that result in spectral lines for some n” 
order nonlinear transformation of the signal. In the case of n = 2, a signal is said to be 
second order cyclostationary if a quadratic transformation produces additive sine-wave 
components that generate spectral lines. This characteristic may be thought of as akin to 
a process being considered wide sense stationary or stationary through order 2. 

To continue the case of n= 2 more specifically, a signal x(t) 1s cyclostationary 
with cycle frequency o if and only if some of its delay products y(t) = x(t)x(t - T) produce 
a spectral line at frequency a. If not all cycle frequencies o for which x(t) exhibits 
cyclostationarity are harmonics of a single fundamental frequency, then x(f) is 
polycyclostationary. Polycyclostationarity implies the existence of more than one 
periodicity in the statistics of a signal. This in turn implies more than one source of 
periodicity such as the case when clocked-digital-data phase modulates a sinusoidal 


carrier in an M-ary phase shift keying scheme [6]. 
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To further illustrate the concept, consider a signal x(t) that contains an additive 


sine-wave component at frequency © and is of the form 
acos(20t +9) with a #0 (3-1) 


The Fourier coefficient defined as 
MY =(x(t)e?™) (3-2) 


where (e) denotes time-averaging, will be non-zero. Note that this is of the form of the 
common representation of the Power Spectral Density (PSD) of x(t) with a spectral line at 
+0. and its mirror at -o. In simple terms, x(f) contains first-order periodicity at frequency 
given by 


IM?| [5(f- a) + 5(f + @)] (3-3) 





Now, considering contributions to x(t) from other sources, the total signal may be 


represented as 


x(t) = acos(2mot +9) +n(f) G4) 


where 7(f) can be thought of as the random energy outside that of the SOI. If n(t) is 
strong in comparison to the SOI such that it masks the sine-wave components in x(t) from 
detection during casual inspection of the waveform, then x(t) can be thought of as 
possessing hidden periodicity. This hidden periodicity can still be exploited by using 
signal processing techniques such as the PSD function as noted above. However more 
powerful cyclostationary signal processing techniques are more sophisticated and can 
unmask periodicities in signals that are hidden even from common traditional techniques 


like the PSD function. 


oA 


B. CYCLIC AUTOCORRELATION FUNCTION 


Traditionally, a common second order, time-domain statistic used in signal 


processing is the autocorrelation function given by the quadratic transformation 
R, (7) = (x()x(t -T)) (3-5) 


In the case of a signal that can be modeled as cyclostationary, this transform will produce 


spectral lines at non-zero frequencies a such that 


Me =(y(t)e?™) #0 (3-6) 


where 
y(t) = x(t)x(t —T) (3-7) 
This signal x(t) as noted above contains second order periodicity manifested in the 


PSD of the delay product given in equation (3-5) above. Transforming (3-5) into a 


symmetric delay product and accommodating complex signals as well gives 


Vio i + < ( -- 4 (3-8) 


which forms the basis for the fundamental second order cyclic moment known as the 


cyclic autocorrelation function: 


a 2s T s Tt —j2na : 
1G {aL =} C =e (3-9) 


Of note, the cyclic autocorrelation function assumes the form of the Fourier 
coefficients of the additive sine-wave components produced by the periodicity of the 


delay product in (3-6). 


Another interpretation of (3-9) is as the traditional stationary autocorrelation 
function multiplied by a kernel e?°™™ that produces spectral] lines at frequencies where the 
stationary autocorrelation function contains additive sine-wave components indicative of 
its periodicity. Consequently, for either interpretation, at a = 0, the cyclic autocorrelation 
function reduces to the traditional autocorrelation function. 

A final interpretation of the cyclic autocorrelation function can be seen by 


defining first 

Wii aije (3-10) 
and 

Vienne (3-11) 


so that u(t) and v(t) are frequency shifted versions of x(t). Now R*(t) can be written 


R-(t)= (: + (i + *)) = R(T) (3-12) 


which is the cross-correlation of the two versions of x(t) shifted up and down by 
frequency OV/2. In other words, the cyclic autocorrelation function may be viewed as the 
correlation in the time domain between two values of x(t) separated in frequency by a. 


Consider as an example, a Binary Phase Shift Keyed (BPSK) signal given by 
ol A.d(t)cos{2nf,t +9) (3-13) 


where d(f) is the binary modulating wave form consisting of positive and negative 


rectangular pulses given by 


2S 


d(t)= ¥d,q(t~t, nf) (3-14) 


ji=-—-oo 


Approximated as a random binary wave, it contains no first order periodicity and thus no 
spectral lines in its PSD. Consequently, the PSD of the BPSK signal s(t) is a scaled sinc 


squared function and will also contain no spectral lines as 1s evident its expression 


ae (3-15) 


Sppsx (f) = an rf T 


Thus, the autocorrelation function of the BPSK signal, the inverse Fourier transform of 


(3-15), will contain no additive sinusoidal components to indicate any periodicity. 


Multiplying the conventional! autocorrelation function by the cyclic kernel e!™ the 
cyclic autocorrelation function follows as 
a I a —j2nat 
Ro =~", (t)cos(2af,T e i 
Bar 
(3-16) 
i ac Sil (Te Platz felt. is mee Cen samen 
b 


where 


Ae) = | a + a Seema (3-17) 


This clearly shows additive sine-wave components at @=+2f.+R,and ~=+k, [9]. 
Thus, the quadratic transform that is the cyclic autocorrelation function unmasks hidden 
periodicity in this simple BPSK signal and proves it polycyclostationary in the process. 
Naturally, the conventional cross-correlation function and the cyclic cross-correlation 


functions are related in a similar manner. 
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: SPECTRAL-CORRELATION FUNCTION 


The frequency domain equivalent for the cyclic autocorrelation function is the 
spectral correlation function (SCF). It follows directly from the Fourier transform of (3- 


9) and has the form 


m l t+Ar/2 : 
Dy a | cp Atlus f FOI2X pu, f - O61 2)du (3-18) 
where 
{+f 6) 
eof Na ie x(uje™!*™ du (3-19) 


is the finite-time Fourier transform of x(t). More commonly used in cyclostationary 
signal processing, the SCF reduces to the conventional PSD at @ = 0 just as the cyclic 
autocorrelation function reduced to the conventional autocorrelation function at @ = 0. 

Continuing the BPSK example from above, the SCF can be found directly from 
(3-17) or by taking the Fourier transform of (3-15). It has the form 


CF) = NG we 
sin= gel lr+s,+£}o (y+ 4,-$] 
a * O —j2m01, 
ro r-5,+}0 (s-1.-%}} 
for a@=n/T, and (3-20) 
a mele a * _ a —J2A(Q+27.\%, 
sneddronBo(r-n-3) 


a ‘ = a —j2a[a-2f,Jt, 
+O f+ f+ Zo y ie se | 


for @=+2f.+n/T, with Q(/) being the Fourier transform of the keying envelope q(t) 


[3]. 


Zs 


The SCF derives its name from the fact that it is in fact a correlation of the SOT in 
the frequency domain. The SCF at frequency fo and cyclic frequency & is merely the 
correlation of two values of the signal in the frequency domain separated in frequency by 


Q) and centered at frequency fo. Figure 3-2 shows this process concept graphically. 


Time Smoothing Method 


, = 
cross-correlation SG) 
yi(t) faerie Te - 


At 


a=, -f; 


Y(v)=X,v+f-a/2), fp=f-a/2 f= Cit 


Vv 


Terminology 


averaging time: At 
cycle frequency: © 
Y,(v) = X,V+f+o/2), fp,=f+a/2 spectral frequency: f 
Vv frequency resolution: Af 
cycle frequency resolution: A a= I/At 
spectral correlation function: S“(f) 





Figure 3-2 - Spectral Correlation Function illustrated, from [10] 


The SCF is plotted on what is called the bi-frequency plane. The plane is defined 
along one axis as spectral frequency f and along the opposing axis as cyclic frequency &. 
Figure 3-3 illustrates this point. The magnitude of the SCF corresponds to a height above 
the bi-frequency plane and is often plotted in that fashion. Note that the values for a 
range from -f,; to f,; while values for spectral frequency f naturally range from -f,/2 to f,/2. 
Because & corresponds to the separation distance of the correlated values in the frequency 
domain, its range extends twice that of f. 

The SCF computation can be highly complex and demanding on any processor. 
For that reason, two algorithms were developed in [7] and designed specifically for 


computation efficiency. The Fast Fourier Transform Accumulation Method (FAM) and 
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Figure 3-4 - SCF of some common modulated signals, adapted from [10] 


yal 


the Strip Spectral Correlation Analyzer (SSCA) both simplify and thereby accelerate the 
computation of the SCF. Because both the SCF and the Cross Spectral Correlation 
Function (CSCF) are necessary for cyclostationary TDOA algorithms, SSCA plays an 
important role in this work. It is presented in more detail in Chapter VI. 

Finally, in addition to their role in the TDOA computations, the SCF and CSCF 
can also be used for purposes ranging from signal detection to characterization and 
estimation of parameters. Specifically, the modulation type of an unknown signal can be 
found using the SCF of that unknown signal and comparing it to SCF(s) of known 
modulation type. Figure 3-4 shows a small library of SCF(s) from several more common 


modulation schemes. 


D. CYCLOSTATIONARY TDOA SIGNAL MODEL 


The signal model used to develop cyclostationary TDOA algorithms 1s very 
similar to that presented in Chapter II. Recall the signals received at two spatially 


separated observers can be given by (2-1) and (2-2) and are below for clarity 


x(t) = s(t) +n,(t) (3-21) 


and 


y(t) = A-s(t— D) +n, (t) (3-22) 


The difference between the two models lies in the temporal and spectral relationships 
between s(t), 2,(¢) and n2(t). In Chapter I, 1,(t) and n2(t) were assumed to have no 
temporally and spectrally components that overlapped the SOI, s(t). In the 
cyclostationary TDOA signal model, )(f) and n2(t) represent all SNOI(s) and noise 
present at the respective observers. They may or may not contain co-channel interferers 
in reality, however for purposes of this work, they are assumed to contain interference 


that in fact spectrally overlaps the SOI. Figure 3-5 portrays this situation graphically. 
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SNOI Observer #1 


Interferer 


Cyclostationary processing provides the 
possibility of eliminating most of the in-band 
interference through the use of cycle frequencies 
unique to the SOI. 


Emitter of Interest 


Interferers reduce SNR and 
SNIR possibly below O dB. 
SNO] 
, Observer #2 


Interferer 





Figure 3-5 - Cyclostationary signal model scenario 


Given this model, the cyclic auto and cross-correlation functions in general can be 


written as 


RY (t) = RE(t) +R, (7) (e223) 
R®(t) =|Al R& (te? + RE (7) (3-24) 
R® (t) = AR®(t—D)e"™ + R®, (7) (3-25) 


just as (2-3) - (2-5) represented the conventional auto and cross-correlation functions. In 


addition, as with (2-6) - (2-8), the spectral correlation functions are 


WES Oem (3-26) 


a 


S*(f)=|A 





Sf emo Sa) (3-27) 


WECM SES Ge SG (3-28) 
Though (3-26) - (3-28) contain the TDOA of the signal D, this parameter is now 

masked by the spectrally overlapping components of m,(f) and n2(t). Thus, any attempt to 
estimate D with traditional methods that compute the above equations at o = 0, like those 
illustrated in Chapter II, would result in a corrupted value. However, if s(t) contains 
some cycle frequency, Mp, that it does not share with any component of 1,(t) and n2(t), 
then cyclostationary techniques can discriminate those contributions to D made by s(t) 
from those contributions of m,(¢) and n2(t) that would otherwise corrupt the estimate. 
Essentially this eliminates the SNOI(s). Thus a reliable estimate of D is possible even in 
highly corrupt environments provided a unique & exists for the SOI. Spectral Coherence 
Alignment is a cyclostationary TDOA algorithm employing measurements of both the 
SCF and CSCF. It determines the TDOA(s) for the TDOA processor developed here and 


appears in detail below. 


E. SPECTRAL COHERENCE ALIGNMENT 


Spectral Coherence Alignment (SPECCOA) was developed and presented in [8] 
using an ad hoc minimum least squares (MLS) approach. Using equations (3-23) - (3- 


25), and the assumption that s(t) contains a unique a, it can be seen that 
Ry (u) = CRY (u- joe (3-29) 


The cross-correlation term for 1,(¢) and n2(t) is eliminated by the use of a cycle frequency 


unique only to the SOI. An estimate of the TDOA, D that minimizes the sum of the 
Squares of error magnitudes between the measured value of the left side of (3-29), RY 
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and the measurement of the right side of (3-29), R® with D substituted yields the MLS 
optimized value for the TDOA. 


Mathematically, from [8], 


pare min} cx()| (3-30) 


Co 


where Ca(T) 1s the estimate of the MLS function. It was further be shown in [8], that the 


solution of the MLS problem is of the form 


Deane man (0) (3-31) 
where 
3 * 
é’(r) =|{ R® (u)R°(u-T) du (3-32) 
=|f[ spss (fy ap (3:33) 


where (3-33) derives from (3-32) through Parseval’s relation. Gardner and Chen go on to 
further prove that (3-33) does indeed peak at t= D. 

SPECCOA derives its name from the fact that the peak in (3-33) occurs by 
maximizing the correlation in f of the two spectral correlation functions through the 
alignment of the phases of the respective functions. 

SPECCOA proves more useful in tactical applications than other cyclic TDOA 
algorithms (see [5] and [9]) by virtue of its ease of implementation and performance in 
corrupt environments [9]. Though other cyclostationary algorithms consistently out 
perform SPECCOA [5], they do so at the expense computational complexity and speed. 
Because this work is intended to demonstrate the feasibility of implementing a 
cyclostationary TDOA algorithm in an operational tactical system, SPECCOA was the 


logical choice for its efficiency and performance. Chapter IV discusses the algorithm 


Sl 


which utilizes TDOA(s) to determine the geolocation of an emitter. Chapter V presents 
the ARL:UT prototype TDOA geolocation system. Chapter VI contains specific 
implementation issues for SPECCOA in the MATLAB® environment. 
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IV. TDOA GEOLOCATION CLOSED FORM SOLUTION 
A. BACKGROUND 


Given the TDOA(s) from pairs of receivers, the problem now lies in determining 
the location of the transmitter from the intersection of the isochrons generated. Many 
distinct methods appear from the geometric interpretation of Schau and Robinson in [1 1] 
to the iterative solution proposed by Loomis in [12] and the closed form solution 
illustrated by Ho and Chen in [13]. While somewhat similar, they all take a specific 
course in determining the location of the transmitter given the determined case. What 
becomes more difficult for each of these solutions is the determination of the contribution 
of errors in the measurement of the TDOA(s) to the final result, the geolocation. 

In [14], Rusu develops a closed form solution avoiding the initial guesswork 
involved in an iterative technique. Unlike [13], his solution 1s completely general. In 
addition, he shows that the errors in the measurements may be propagated through the 
mathematical model in a linear fashion, simplifying the determination of uncertainty in 
the geolocation. The following development draws exclusively from his derivation. 

The problem of determining the location of an emitter in three dimensions given 
the times of arrival (TOA) at four spatially separated receivers has often been treated as a 
TDOA problem. That is, time is treated as absolute time. However, as Rusu points outs, 
the mathematical model is invariant to time translation and thus may be treated as a TOA 
problem by setting the arrival time at any one of the receivers as the origin of the 
problem, t = 0. This simplification without loss of generality facilitates the development 
of linear error propagation from the initial measurements through the mathematical model 
to the final result as noted above. 

In the presence of moving emitters or receivers, both TDOA and FDOA must be 
determined by one process as shown in the development of the CAF in Chapter II. The 
case of FDOA can be similarly simplified by assigning a Doppler shift of zero to one of 
the receivers reducing the problem to one of frequencies of arrival (FOA). Fortunately, 
given the independence of the TOA portion of the model from the frequency of 


transmission and Doppler shifts, the TOA and FOA equations may be treated separately. 
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Though TDOA is the emphasis of this thesis, for completeness both the TOA and the 


FOA solutions will be developed here. 


B. MATHEMATICAL MODEL 


The explicit functional model relating the observable X to the dependent variable 
of interest Y, in the TOA case is the solution of the set of scalar equations of the form 
F(X,Y) = 0. F may be referred to as the implicit functional model, X the independent 
variable and Y the dependent variable. In the development of the TOA-FOA problem, X 
is a 32-dimensional real vector state space of observations while Y is an 8-dimensional 


Status vector state space. Figure 4-1 illustrates the scenario in general. 


receiver | 
heecilver 2 State vector: (X), Yp Zp Mp Yp Wy fy, Oy) 
State Vector: (Xz, V2, Za, Uy Vay Wo, by, Wy) 


target 
State vector: (x, y, Z, u, Vv, W, t, @) 


receiver 4 


receiver 3 
State vector: (Xp V4 Zp Uy Vs, Wa ly Wy) 


State vector: (3, Vz Zy Uy Va, W3, fy, @) 





Figure 4-1 - TDOA model 


The TOA-FOA case assumes that the i-th receiver located at rj = (Xi, yi, Zi) and 
moving with velocity v; = (uj, Vj, W;) detects a signal at time ¢; with frequency @,. Four 


observers for the determined 3-dimensional case, constitute the observation state space 
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vector of 32-dimensions. The parameters of the transmitter, its location and velocity 
given by r = (x, y, z) and v = (u, v, w) respectively, its transmit time ¢ and frequency @ 
compose the 8-dimesional status vector state space. 

The fundamental equation relating the signal travel time to the separation distance 


between observer / and an emitter assuming the signal travels at the speed of light is 
x, —2) +, +2) =e, -9 (4-1) 


where 4/(x; — x) + GQ — y)? Pia z)’ is the Euclidean distance between r, and r denoted 


hereafter by 








= r|| and c is the speed of light. Forming appropriate vectors using a 


single equation to include all four observers produces the implicit functional TOA model 
F(X, Y) =r, -r||— c(, - 1), aie (4-2) 


Simple algebraic transformation leads to 


(x, - x)? +(y, -y)? +z, - 2)” -(,-1)’ P =0, ee ee) 


In vector form the TOA equations are defined as Fy = (F), Fo, F3, F4). 
The Doppler equations can be similarly treated. The fundamental relationship 
between the Doppler shift observed by receiver n and the relative velocities of the emitter 


and observer is 


a -of =*) an (4-4) 


c) Ir, —4 








where @4q is the Doppler shift observed, @ is the transmitted frequency, (Vv, - Vv) 1s the 


Cae 
relative velocity and _— is the unit vector from the observer to the transmitter. 


r 
r, —r| 
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Again forming vectors based upon this relationship yields the implicit FOA functional 


model 





F,(X,Y)=@, -of 1 =): = i = 1,2,3,4. (4-5) 


defined as Fp = (Fs, F6, F7, Fs). 
Combining the TOA and FOA models produces the TOA-FOA implicit functional 
model defined as F = (Fy, Fp). 


ve, 


C. DETERMINATION OF THE TRANSMITTER STATE VECTOR 


The determination of the 8-dimensional transmitter state vector involves finding 
for X, the parameter Y that satisfies the combined TOA-FOA implicit functional model, 
F,(X,Y) = 0, i = 1,2...,8, where the functions F; are defined in (4-3) and (4-5). The 
solution of such a system of equations occurs in two steps. First, the TOA equations can 
be solved for Yy= (r, f). The resulting solution can then be used in the FOA equations to 
solve for the unknowns Yr= (Vv, @). 

The TOA equations are irrational in their form in (4-3) and thus must be squared 
in order to be solved. This produces a set of quadratics for which two sets of solutions 
exist corresponding to the two roots of the equations. In most cases, Rusu points out, the 
two solutions to the rational quadratics are also solutions to the original irrational set of 
equations. Rarely, most often in the 2-dimensional case, only one of the two solutions 
also solves the irrational set of equations and leads to a unique solution to the geolocation 
problem. 

Since each solution to the TOA case also produces a unique solution to the FOA 
equations, more often than not, there are two distinct solutions to the problem. This ts the 
classic case of ambiguity encountered with every TDOA solution proposed thus far. 

Rusu handles this problem by noting that information outside the algorithm must resolve 
the ambiguity. Prior knowledge of the general area of the target’s position might 


eliminate one of the solutions. Multiple measurements may also reveal one solution 
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converging in one location while the other solution diverges and produces seemingly 
unrelated geolocations. 


To solve the TOA equations, (4-3) must be squared to produce 


(x, -x)°+(y,-y)° +(z,-z)° -(t, -1)° =9, i = 1,2.3,4. (4-6) 


The solution to (4-6) is found by subtracting the equations two-by-two to eliminate x’, y’, 


z and t?. This produces a set of linear equations in r and ¢ that may be written as 
Ar=qt+s Ca 


As the equations can be subtracted in any order to produce such an elimination, 
there exists many possible resulting systems of equations. One more obvious possibility 


used by Rusu is 


X,7~%X_ VY, 7 V2 <p 7 


A=|x,-X; Y,-y¥; 2-2; (4-8) 


X37~%q V37 Vg <3 7 hs 


th, 
f= |e ie ae) 
ety 


(hf -faf-2 +2 
S= 5 lI = lrI| —t; +1; (4-10) 


Ps 2 2 5 
Isl - hall - 3 +4 


This system makes it possible to solve for r in terms of t producing a f- 


parametrized solution to the TOA equations 


r=gt+h, (4-11) 


So 


where 


g=A'q and h=A's. (4-12) 


The introduction of this result into the k™ range equation from (4-6) yields a 


quadratic in time t. This equation, with some algebraic manipulation can be written as 


at’ + 2bt+c=0 (4-13) 
where, 
a=|el -1 (4-14) 
b=g-h+g-r, +1, (4-15) 
c=|h-r,° -2 (4-16) 


The two roots that are the solution to (4-13) produce two values for r when 
inserted into (4-11). Rusu notes that the roots are usually real but in the rare cases when 
they are complex, experience has shown him that the real part of the complex roots 
suffices as a solution. 

Once the TOA equations provide the emitter position and time of transmission, 
these values produce the transmitter’s velocity and transmitting frequency from the FOA 


equations. If the FOA equations are rewritten as 


l 
Sy aie one i= 1,2,3,4, (4-17) 
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the FOA case becomes a non-homogeneous linear system that may be solved using 


standards procedures like Gauss-Jordan elimination [14]. 


D. ERROR PROPAGATION 


Perhaps the most valuable portion of Rusu’s solution comes with the error 
propagation. He developed a linear method through which errors in the measurement can 
be propagated through the model and result in predictable errors in the fix. The following 
derivation is based solely on his development of this method. Though he develops both 
the TOA and the FOA propagations, only the TOA case will be presented here. 

The key to the linear propagation of errors is based upon the assumptions. First, 
Y = Y(X) is assumed to be differentiable in the neighborhood of the measured value X 
containing the exact (error free) value of the observable X. Second, if the errors in the 
observation and the parameter are be denoted dX = X - X and oY = Y - Y, Y is assumed 


to be linear in the region of these errors and thus oY can be written 


OY = (Spx (4-18) 
OX 


The following Jacobian matrix represents the derivative of a vector valued 


function, Y = (¥1,Y9,...,%,) with respect to a vector variable X = (%1,%9,....Xn) 


le or 4-19) 


A1Vp X2Vn Xm\n 


Given the linear approximation in (4-18), the covariance matrix of the parameter 


can be related to the covariance matrix of the observation through 


62 


dY oY 
c -[™)o(™ 20 


The implicit function theorem provides a convenient method for computing the 
derivatives of Y with respect to X from the relation F(X,Y) = 0. This provides the 


general formula 


aay “ 
OX OY) \dX 


Just as the relationship of the TOA and FOA solutions allowed the separation of 


4D 


the TOA case from the FOA case, the matrices involved in (4-21) can be written such that 


the TOA and FOA contributions are separated as 








Mp ON 
ON Slag Xemerc Xe 
OX | A%e OY, “a 
dX, aX, 
oF ak 
One ox 
OX | aR aR —<— 
OX, aX, 
yg 
OF _| oY, (4-24) 
OY | Fp OFp 
aY, oY, 


The independence of the TOA solution from the FOA equations is evident in the 


. OF OF 
zeros found above corresponding to x. and why ; 
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The problem of inverting the Jacobian, dF/dY, can be solved numerically 
however, Rusu managed to express the inverse of dF/0Y in terms of the inverses of the 
TOA-FOA blocks forming it. He points out that this formally separates the TOA and 


FOA error propagations and allows computation of just the TOA solution and errors if 


just the transmitter location is desired. The block inverse of (4-24) is 


OF, ) : 
by . Mr =) =i (4-25) 
oY (Fe a | Ea 


oY, } OY, \ oY, oY, 











It is now possible to write the components of the right hand side of (4-22) as 



































=| 
OY, OF, OF, 
=— 4-26 
OX, eo SS ) | 
OY, _ 0 (4-27) 
OX, 
Ove | Ob, (dF, | aF, ) { oY, _{ OFe : OFy (4-28) 
OX, 7 OY, ane dY, OX, OY; OX, 
I 
OYr _{ Fp ) ( OFe (4-29) 
ONpeeee ype) Co 


using (4-21), (4-23) and (4-25). 
Now the TOA Jacobian matrices related to the TOA error propagation can be 
explicitly expressed using the functional model given in (4-3). The Jacobians can be 


expressed as 


4] 


Xo, Yor 2% Zoi 


oF, Xoo Yor 22 Lop 


= ' (4-30) 
OX, Xo; Yo3 S23 SL Q3 
Xos Vos %og Log 
Xo Vine cio 0 0” One 0 O O 
OF, : 0 0 0 “O° Xe ysn cope ea 0 Onea0 (4-31) 
ab 0 0 O° OALEO 0 OG = 0 0 O O 


X39 Yao 240 fog 


where 


xij = Mes 


Implementation of Dr. Rusu’s algorithm into MATLAB® followed naturally 
given its matrix and vector nature. It plays an integral part in the ARL:UT prototype 
TDOA geolocations system which 1s discussed in Chapter V. Specifics regarding the 


encoding of Dr. Rusu’s algorithm can be found in Chapter VI. 
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Nie THE CARRY-ON MULTI-PLATFORM GPS-ASSISTED TIME 
DIFFERENCE OF ARRIVAL SYSTEM 


A. PROJECT BASIS 


Given the problem of developing a portable, GPS assisted TDOA system with 
commercially available products, the project team at ARL:UT set out to define the 
problem and develop solutions based on which portions of the problem could be most 
influenced. The descriptions and performance reports below are based upon [15] and 
numerous conversations with team members. 

The TDOA problem can be reduced to three sub-problems: geometry, signal 


timing and measurement. All three are related to the rms error of the geolocation by 


= GDOP:-o 


O position measurement iS =| ) 

The Geometric Dilution of Precision (GDOP) is a scalar multiplicative factor 
which serves to magnify any timing and measurement error. It is directly related to the 
relative positions of the observers and the emitter. Poor GDOP(s), generally numbers 
greater than three in this application, result from geometries where the observers are lined 
up or clustered in one area with respect to the emitter. Conversely, good GDOP(s), 
numbers smaller than about three, result from situations where the observers are 
positioned such that their TDOA(s) intersect at nearly right angles. In three dimensions, 
the most favorable GDOP occurs by maximizing the volume of the polyhedron formed by 
pointing four vectors from the four observers to the emitter and enclosing the figure that 
results. 

The geometry of the situation is essentially out of the control of the TDOA system 
as the location of the emitter is not known a priori but would be needed in order to 
position the observers in an optimal configuration. Only when many more observers 
participate than the minimum needed for a determined solution (three are needed for 2D 
and four for 3D) can the geometry of the situation be better controlled by the TDOA 
system. In practice, this occurs rarely and thus, the best course of action considering (5-1) 
above is to minimize the timing and measurement error to the greatest extent possible. 
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GPS provides the capability to reduce both these error sources significantly. It was with 
this in mind that the ARL:UT project team, led by Mark Leach, pursued the design and 
implementation of a new GPS-assisited TDOA system. 

NSGSA in its tasking to ARL:UT added some additional requirements in the area 
of component sources, architectures and interfaces with existing equipment. The 
following sections outline the system in progressively more detail starting with a general 
summary and concept of operations, continuing into a hardware description and ending 
with a more detailed discussion of error sources and system performance in operational 


tests. 


B. PROJECT OVERVIEW 


v 


The system would receive only a frequency of interest in the HF, VHF or UHF 
line of sight communications frequency bands. It must utilize an open system 
architecture and maximize the use of commercial and government of the shelf hardware 
(COTS and GOTS respectively). It was permitted just a single voice grade channel for 
communication between observers and needed to interface with existing Navy hardware 
for that purpose (i.e. AN/WSC-3 transceiver). Finally, the system needed to be 
compatible with the Navy’s Unified Build (UB) environment and interact seamlessly with 
the Joint Maritime Command and Information System (JMCIS). 

The system in its current configuration consists of a network of observers 
(minimum of three for a geolocation in 2D and four for 3D) with one observer acting as 
the master node, the others as slave nodes. It employs the distributed processing 
technique depicted in Figure 5-1. Once given the frequency of interest by an existing 
signal acquisition system , the master node tunes its receiver to that frequency and 
notifies the slaves to do likewise. Each node samples the incoming signal and buffers the 
data in mass memory storage. The master logically determines the 400 ms portion of its 
samples that contains the most energy from the SOI and computes the Fast Fourier 
Transform (FFT) of this 400ms. The Master then broadcasts the FFT to the slaves with a 
GPS time stamp. 
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Figure 5-1 - Distributed processing, from [15] 


Once the slaves receive the FFT from the master they search their buffers for the 
400 ms of samples that correspond to the epoch of the time stamped FFT. They then take 
the FFT of their data and perform a CAF to extract the TDOA. In addition to this TDOA, 
the slaves send their GPS derived position, time and velocity (PVT) measurements 
corresponding to the time of intercept to the master. 

The master collects the TDOA(s) and PVT(s) and forwards the package to a 


standard Navy TAC3/4 workstation which calculates the geolocation. 


C. HARDWARE DESCRIPTION 


The system 1s comprised of two major sub-systems: the measurement subsystem 
and the geolocation subsystem. Each is designed in modular fashion accepting certain 
formatted inputs and producing certain formatted outputs. This modularity facilitates 
changes internal to either system including entire components or software modules as the 
process within each of these subsystems is inconsequential to the other module. Only the 
inputs and outputs are relevant. 

The measurement subsystem includes the hardware and software used to compute 
the TDOA(s). It is based on a 13 slot, C-size VXI chassis housing an HP Signal Analyzer 
package. This package includes a Motorola 68040 based VXI controller, a 10 MHz, 23- 
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bit (18 linear) analog to digital converter (A/D), a digital signal processing (DSP) card, 
200 MB of disk storage and software. The DSP module is a Motorola 96002 based 
floating point DSP chipset on a VXI card. Also on the chassis are a Ball-Efratom model 
FRK 10 MHz Rubidium (Rb) reference oscillator, a Precise Positioning Service (PPS) 
capable GOTS GPS receiver by Ashtech Inc. and E-systems and an ARL:UT developed 
clock synchronization module. The phase lock loop (PPL) module allows the Rb 
reference oscillator to be in synch with the GPS | pulse per second (pps) timing output. 


The bulk of this subsystem is shown in Figure 5-2. 
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Figure 5-2 - Hardware configuration, from [15] 


The geolocation subsystem includes a TAC3/4 computer on which reside the 
algorithms both to compute the geolocation and communicate that fix to the JMCIS. In 
addition, it also holds the system’s capability for simulation and playback modes. The 
TDOA algorithms consist of the closed form solution presented in Chapter V and a 
Kalman filter method. The closed form or determined solution is the primary solution 
and lends itself quite naturally to situations where the emitter has short transmission 


durations which preclude the use of the Kalman filter method. In the event of long 
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duration transmission, the determined solution will feed the Kalman filter a first 
prediction and allow the Kalman filter to produce fixes as a function of time over the 
duration of the transmission. The Kalman filter, when given sufficient data and a starting 
point, has better error prediction capabilities and is thus the preferred method in the 
correct situations. 

The interface with JMCIS provides the system a data display capability. By 
updating the track database through a standard UB message, the system is able to both 
communicate with the Navy’s latest command and information system and avoid having 
to produce its own man-machine interface on this output level. 

Finally, this system.is also capable of simulating scenarios in order to investigate 
the effects of variations on the many factors involved from the receipt of the signals to the 
processing of the data to the computation and display of the fix. In addition, a capability 
to playback sessions recorded and stored on tape 1s built in to the geolocation subsystem. 
This feature allows for post processing and analysis of operational tests enabling finer 
detailed analysis of many of the variables including error source investigation. 

Perhaps the one key point that stands above the other accomplishments of the 
system is the significant reduction of signal timing errors between the collection nodes. 
With the use of the PPS-capable GPS receivers, the 1 pps output from the GPS satellites 
in concert with differential GPS techniques enables the time offset between two GPS 
receivers, separated by up to hundreds of kilometers, to be on the order of 20 ns. In 
addition, the PPL circuitry allows each node’s Rb standard to be synchronized to within 5 
ns of this 1 pps output giving a total signal timing error of a mere 25 ns. In terms of light 
distance, that translates to 7.5 meters. Figure 5-3 shows the entire system’s functional 


relationships in summary. 
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Figure 5-3 - Tactical TDOA Functional Diagram, courtesy ARL:UT 


D. ERROR ANALYSIS 


Many sources of errors exist in the operational system. Another accomplishment 
of the project is the detailed identification and prediction of the magnitude of those error 
sources which have been validated by the operational tests outlined in the next section. In 
their investigations, the ARL:UT team was able to collect the sources into four major 
categories: fundamental limits, instrumentation limits, operational environment and 
network geometry. 

The lower bound on timing and measurement accuracy may be related to the 
variance of the TDOA estimator. This in turn is related to the incoming SNR, signal and 
noise bandwidth and integration time as first shown by Hertz and Azaria in [2] and later 
extended more specifically to this system by Rusu and Giulianelli in [16]. 

Instrumentation limits add significantly to the errors in signal arrival time 


determination. Specifically, though ideally each GPS receiver will output the GPS 1 pps 
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at the same time during each one second epoch, there will be some offset in reality. In 
addition, that 1 pps timing pulse is used to trigger the sampling of the incoming signal. 
Delays from the time the A/D converter receives the 1 pps and the onset of sampling for 
that clock cycle will vary from one A/D to another. This compounds the | pps offset 
error and translates directly to error in the cross-correlation computation of the TDOA(s). 
The cross-correlation methods of TDOA determination measure the time difference by 
determining a peak in the computed function. Five unknowns exists in this computation: 
signal frequency, signal bandwidth, noise bandwidth, SNR, and statistical properties of 
both signal and noise. The width of this peak and the presence of multiple peaks in 
indicative of the magnitude of the lower bound on the error in the measurement discussed 
above and 1s the direct result of a combination of these factors. 

The operational environment and network geometry constituted the bulk of the 
uncontrollable error contributions. While during testing, the environment is somewhat 
controlled, the operational system would certainly be subject to the conditions of the 


current tactical situation. This holds true as well for the geometries. 


the PRELIMINARY TEST RESULTS 


Published test results include both static and dynamic emitters as targets. In 
November 1994, the geolocation of a static emitter in 2 dimensions was demonstrated in 
the Austin, Texas area using three observers with baselines of approximately 33 km. The 
target was a Stationary NOAA FM transmitter at 162.4 MHz. The system at this point 
tested with a 2-D rms error of 124 m. Details can be seen in [15]. 

Dynamic tests were conducted in January 1995 using a moving emitter, two 
stationary observers and one moving observer. Again the tests were conducted in the 
Austin area and in just two dimensions only given just three observers. The target emitter 
was an aircraft at approximately 2000 feet, traveling at speeds of approximately 100 knots 
and transmitting with a 20 W UHF signal at 461.1 MHz. The moving observer was a van 
mounted node traveling at speeds of up to 35 mph in a 2 square mile area. During this 


test, the prototype system performed geolocations to within 277 m rms error using the 
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Kalman filter solution. This accuracy was achieved despite various geometries that 
ranged from outstanding GDOP(s) of almost | to very poor GDOP(s) of 10. 

While these tests failed to achieve the goal of 100 m rms error, instrumentation 
changes promised to increase the accuracy without significant changes to the system. 
Two additional tests, each more extensive and demanding than the previous, have been 
conducted with the system. While no performance results have been formally published, 
shipboard tests off the coast of San Diego California have proven the system capable of 
120 mrms error against representative VHF/UHF dynamic emitters with excellent 
geometry and ~400 mrms error with poorer geometries [17]. Further changes in 
receivers promise to remove biases in the measurements created by the current radios and 
improve performance proportionally. 

In addition, the use of a cyclostationary TDOA determination algorithm should 
provide some improvement in the performance of the system low SNR environments. 
The initial step for the incorporation of cyclostationary signal processing into the 
ARL:UT system is to code the algorithms in MATLAB® and test them with ARL:UT test 
data. The coding of the algorithms in MATLAB® appears 1n Chapter VI. 
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VI. ALGORITHMS 


A. BACKGROUND 


Any implementation of a TDOA processor utilizing SPECCOA must provide both 
an SCF as well as a CSCF at acycle frequency 0, characteristic of the SOI. In [19], 
Loomis and Berstein develop a functional TDOA processor using SPECCOA that 
provides for the efficient computation of both spectral correlation functions as well as the 
selection &. This system represents the basis for the implementation of SPECCOA and 


the geolocation algorithm in MATLAB® for this project and appears in Figure 6-1. 
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Figure 6-1 - TDOA geolocation processor functional diagram, adapted from [19] 


The four blocks of the TDOA processor constitute progressively lesser degrees of 
computational complexity. The most complex of the operations is the calculation of the 
SCF(s) in the first operation of the processor. The SCF(s) of the master and the four 
slaves are necessary for the selection of the appropriate cycle frequency in the next 


operation of the processor, the characterization. The SSCA algorithm computes the four 


i 


SCF(s) and uses a plotting routine to display them. In the next operation, the 
characterization 1s performed manually using a plot of the maximnum values of the SCF(s) 
for each value of cycle frequency. This plot facilitates the choice of the optimum cycle 
frequency for use in SPECCOA. The third operation involves the computation of the 
SCF and CSCF at the cycle frequency selected during characterization. Computed with a 
simple frequency smoothing algorithm, the SCF and CSCF at the cycle frequency of 
interest feed the final TDOA processor operation, SPECCOA. SPECCOA in turn 
determines the TDOA by plotting the solution to the ad hoc MLS optimization found in 
(3-32). The maximum value of this function represents the TDOA. 

The SCF(s) in the first operation as well as the characterization must be computed 
for all four observers for each set of data. Cycle frequencies for the same emitter vary 
with observer due to instrumentation mismatches at each location. Details of cycle 
frequency selection in these situations appear in Chapters VII and VIII. The cycle 
frequency specific SCF and CSCF are computed just three times because as noted in 
Chapter V, the only TDOA(s) used 1n practice are those between the master node and the 
three slaves. No cross slave TDOA(s) are computed or used in the prototype system and 
thus they are not computed here. Finally, three runs of SPECCOA are required to feed 


the geolocation algorithm. 


B. STRIP SPECTRAL CORRELATION ANALYZER IN MATLAB® 


The SCF(s) computed in the first operation of the TDOA processor can place a 
significant burden on the processor 1f not optimized for computational efficiency. In 
order to achieve the efficiency necessary for use in a tactical environment, a highly 
expeditious method for computing and plotting the SCF is used. Directly from its 
definition, (3-17), the SCF can be computed using a traditional time-smoothing approach 
as depicted in Figure 6-2. While this approach provides the most accurate estimate of the 
SCF, it proves prohibitively time consuming and inappropriate for most applications [20]. 
The SSCA algorithm, a modification of the time-smoothing method, provides the best 
combination of accurate estimation of cycle frequencies of interest and computational 
efficiency [21]. By eliminating one of the band-pass filters in Figure 6-2 and using a 
Fourier transformer at the input and output, SSCA computes the SCF along 
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Figure 6-2 - Time-smoothing spectral correlation analyzer [20] 


diagonal strips of the bi-frequency plane as depicted in Figure 6-3. While the SSCA does 
manage to compute the SCF in the most efficient manner developed thus far, it has some 
compromise. The output signal to noise ratio suffers slightly as a result of the efficiency 
gained by the elimination of the filter and the use of the Fourier transform. However, the 
computational savings gained far outweigh the minor degradation in signal power to 
noise power at the output of the analyzer. 

Figure 6-4 illustrates the SSCA architecture. The input filter consists of a sliding, 
Hamming-windowed coarse FFT of length NV’. Each of these N’ length segments is 
essentially the output of a band-pass filter with a bandwidth of 1/ N’. This band-pass 
output is then downconverted to baseband by a complex exponential multiplication. 
These products are termed the complex demodulates of the input signal and appear in 


greater 
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Figure 6-4 - SSCA architecture, from [20] 
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detail in Figure 6-5. As can be seen, the input sequence is essentially decimated by the 
process of windowing and sliding at the input filter. Loomis and Brown found that the 


most appropriate value for the decimation factor or subsampling parameter L was N’/4 


[20]. 
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Figure 6-5 - Detail of the computation of complex demodulates 


Because these complex demodulates are multiplied by the original sequence 
before the final FFT, they must be interpolated to the original sampling rate. In [20] and 
[21] Loomis and Brown use a simple hold operation. Linear interpolation induces less 
high frequency error and adds little in computational complexity. It 1s used here in place 
of the first order sample and hold interpolation. Finally, the interpolated complex 
demodulates and the original sequence are complex multiplied and the product Fourier 
transformed by a full length, unwindowed FFT. The relationship of the output to the bi- 
frequency plane appears in Figure 6-6. Each row of output represents one of 


the N’ diagonal strips in the bi-frequency plane as depicted on the previous page in Figure 
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6-3. Each column represents the values of the SCF along the diagonal lines. The 
resolution along lines of fis coarse at 1/ N’ as a byproduct of the computational savings. 
However, the resolution in a, at 1/N, is much higher which accommodates the 


characterization nicely. 
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Figure 6-6 - Block form of SSCA output 


Given that the majority of operations in the SSCA can be implemented as vector 
Or matrix operations, its implementation in MATLAB® 1s fairly simple. The psuedocode 
for the algorithm appears in Figure 6-7. The MATLAB® code for the computation of the 
SCA over the entire bi-frequency plane, called ssca, appears in Appendix A. Of note, 
several MATLAB® functions forced the use of additional statements that might not 
otherwise be found given the psuedocode. The fft function in MATLAB® computes both 
positive and negative frequencies but does not place the origin in the center of the output 
vector. Instead, the output of fft takes advantage of the cyclic nature of the DFT. Its 
output begins at the origin of the frequency, moves through f,/2 and then mirrors the 


negative frequencies with the positive continuing at -f,/2 and ending with zero again. 
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This forces any application needing a double sided spectrum to rearrange the vector so 
that the negative frequencies start the output rather than end it. MATLAB® has an 
fftshift function that performs such an operation and appears in the code in several places 
for that reason. The output of ssca is a matrix with N’ rows corresponding to 

the N’ diagonal strips of the computation and columns corresponding to the SCF along 
the respective strips. Referring again to Figures 6-3 and 6-5, these strips are not 
computed along lines of f and o and must therefore be transposed along those lines before. 


plotting. 


/* Compute Complex Demodulates of x */ 
Do p:= 0 to P-1 
Compute x,(pL,f,) = FFT[a(r)x(pL+r)] 
Dok := —N’/2 to N’/2-1 
Compute X,(pL, f,) = xT(pL, f, )exp{-j2mpLk/N  } 
end 
end 
/* Interpolate X,(pL, f,) */ 


Compute X,(pL, f,) = X,(n, f,) 


/* Compute and Smooth Product Waveforms 
Dok := -N’/2 to N’/2-1 
Compute Si fh) = CAG 


Compute S% (n, f,),, = FFT{g(n)S/ (n, f,)} 
end 





Figure 6-7 - Psuedocode for SSCA computation, from [21] 


The program plotscf performs that function transposing the matrix output of the 
ssca program and plotting both the 3-D SCF and the 2-D magnitude of SCF versus «& for 
the characterization operation. The transformation involves applying the relationship 
between the diagonal strip subscript k, the frequency f and the cyclic frequency a. Those 


relationships from [22] are 


2 











5-4) sks a1 (6-1) 
a=2f,-2f (6-2) 


Given these relations and the MATLAB® surf function, the SCF plots easily 
along lines of fand a. The 2-D plot for additional characterization aid follows directly 


from the SCF plot. 


C. SPECTRAL COHERENCE ALIGNMENT IN MATLAB® 


SPECCOA represents yet another algorithm highly suited for MATLAB®. The 
majority of functions associated with computing (3-32) are vector operations which 
simplifies the MATLAB code considerably. Figure 6-8 1s a basic block diagram of 
SPECCOA. 
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Figure 6-8 - SPECCOA block diagram 
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Once the characterization has determined the cycle frequency of interest, Ot, a 
frequency smoothing algorithm estimates the SCF and CSCF for G&. The inverse Fourier 
transform of the complex product is then multiplied by the kernel e’”*”. The maximum 
value of the real portion of this product represents the TDOA of the SOI between the 
observers. 

SPECCOA is used three times to determine the TDOA between the master node 
and the three slave nodes. These TDOA(s) along with observer positions are inputs to 


Dr. Rusu’s geolocation algorithm the output of which is the location of the emitter. 


D. THE TDOA CLOSED FORM SOLUTION IN MATLAB® 


Dr. Rusu’s closed form solution lends itself readily to MATLAB implementation. 
The inputs include four observer position vectors in WGS84 coordinates and three 
TDOA(s), D2, Di3 and Dy4. All the intermediate variables are either vectors of length 3 
or matrices no larger than 3 x 3. The final solution involves finding the roots of a 
quadratic that represent the time of emission of the SOI and substituting that time into a 
range equation to find the WGS84 coordinates of the transmitter. 

Of particular note to the MATLAB® code that appears in Appendix A is the use 
of the backslash function, \, in lieu of the matrix inverse in the algorithm. This function 


was specifically developed for use in solving equations of the form 


aed (G55) 


for g where A is a matrix and g and q vectors. Normally, in matrix algebra the solution 


has the form 


g=Aq (6-4) 


if A has an inverse. However in MATLAB®, it is more accurate to use the matrix 


division function or backslash. The accuracy of the backslash in this scenario approaches 


Se 


machine precision while the use of the matrix inverse can reduce that accuracy several 
orders of magnitude [23]. 

One problem persists with the algorithm but is a function of the closed form 
solution rather than its implementation. The two roots both represent viable solutions to 
the problem. Information outside the algorithm must be used to eliminate the ambiguity. 
Multiple solutions from moving observers would more often than not cause one solution 
to converge while the other diverged. No matter how its resolved outside the algorithm, 
the ambiguity always poses a problem insurmountable by the mathematics. 

The algorithms need to be tested on controlled data prior to evaluation with 
ARL:UT test data. Several signals were generated using Simulink®. These signals serve 
to ensure the proper functioning of all the code written for this work. Following the 
generated data, ARL:UT data provides a vehicle for comparing the current CAF based 
ARL:UT method with SPECCOA. The test plan is detailed in Chapter VII. 
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Vil. TEST PLAN 


A. GENERATED SIGNALS 


Three signals are used in initially testing the TDOA processor. Generated with 
Simulink®, all three are BPSK signals and vary from noiseless to very poor SNR. Each 
signal was sampled at 100 kHz. They had carrier frequencies of 0.25 f, and data rates of 
0.05 fs. All were modulated with a random binary signal generator. They are used in 
various combinations to verify the functions of the three major programs, ssca, plotsef 
and speccoa. 

The first set of tests involves using a noiseless BPSK signal to first verify ssca 
and plotscf. Given that these programs function properly, speccoa is tested finally. 
Figure 7-1 shows the simulation used to generate noiseless BPSK. A total of 1 second of 


signal is actually collected and stored as testbpsk.mat. 


Random Binary 


Wave BPSK Modulation testbpsk.mat 





Figure 7-1] - Generation of noiseless BPSK as testbpsk.mat 
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This signal approaches the ideal for the processor for several reasons. First, it is 
completely binary from its inception. No analog to digital conversion is necessary and 
thus introduces no error in the process. Second, it is completely noiseless. Finally, it is 
modulated by a random binary wave which introduces no biases in the modulated 
waveform. Given this, its SCF should approach that given by theory taking into account 
the time smoothing effects of the SSCA algorithm. Figures 7-2 and 7-3 are the time 


domain and frequency domain plots respectively for testbpsk.mat. 
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Figure 7-2 - Time domain plot of testbpsk.mat 


The ideal SCF of a noiseless BPSK signal can be seen in Chapter III, Figure 3-4. 
The data feature of testbpsk.mat should appear at 0.05 f; while the twice carrier feature 
in conjunction with the data feature appear at 0.45 f,, 0.5 f6 and 0.55 f,. Figure 7-4 
displays the results of ssca and plotsef for N = 4096 point sample. Only the positive 


cycle frequency portion of the bi-frequency plane is plotted. The cycle frequency axis 
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draws from upper right to lower left and is plotted in terms of nf,/N, n as axis label, I/N as 
cycle frequency resolution. The spectral frequency follows from upper left to lower right 
labeled in the same manner. They are unlabeled in the figure to eliminate the additional 
clutter. As noted in Chapter VI, the SSCA algorithm trades resolution in the frequency 


domain for computational efficiency without losing the four key features. 
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Figure 7-3 - Frequency domain plot for testbpsk.mat 


Finally for testbpsk.mat, speccoa is run with observer | collecting testbpsk.mat 
samples 1 to 8192 and observer 2 collecting testbpsk.mat samples 71 to 8263 for a 70 
sample delay. The results of this test appear in Figure 7-5. Because the signals at the two 
observers are completely correlated, this scenario represents an ideal situation for the 
algorithm. Its complete success is expected for this case but serves as an initial check of 
speccoa. The remaining two signals introduce both noise and poorly correlated signals to 


the scenario. 
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Figure 7-4 - SCF (top) and Mag(SCF) vs alpha (bottom) for testbpsk.mat 


computed and plotted by sscea and plotscf 
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Two very noisy BPSK signals constitute the remaining generated signal testing for 
the system. They were created by passing the noiseless BPSK signal through separate 
AWGN channels with E,/No of 0.88 or -0.534 dB. This corresponds to probability of bit 
error of 107. Figure 7-6 shows the generation of nl bpsk.mat and n2bpsk.mat. Because 
the AWGN channel introduces significant noise the BPSK pena nlbpsk.mat and 
n2bpsk.mat are highly corrupted. Figures 7-7 and 7-8 show the time domain and 


frequency domain plots respectively for the two signals. 
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Figure 7-5 - Output of speccoa routine for testbpsk.mat with 70 sample delay 


Naturally, the level of noise in these BPSK signals will be evident in their SCF(s). 
Figures 7-9 and 7-10 display the output of ssca and plotscf. Note that the reduction in 
SNR compared to testbpsk has added noise to the plots but has not masked the important 


features. While in theory as presented in Chapter III, the AWGN should have no features 
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other than those evident at © = 0, it does add noise to the higher values of ain the SCF 
plots in Figure 7-9. It does so because the AWGN channel used in the simulation is not 
ideal AWGN and thus is not truly random and Gaussian. Thus, some statistical 
impurities in the channel as well as the computational trade-offs in the estimation 


procedure of the SSCA algorithm lead to contributions above a = 0. 
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Figure 7-6 - Generation of corrupted BPSK as nlbpsk.mat and n2bpsk.mat 


Finally, speccoa uses nlbpsk.mat for observer | and n2bpsk.mat delayed by 70 
samples for observer 2 to verify that it indeed works in highly corrupt environments given 
the proper signal features as inputs. The results appear in Figure 7-11. The successful 
conclusion of testing with generated signals leads next to testing of the ARL:UT data 


from a November 1995 test in the Austin, Texas area. 
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Figure 7-8 - Frequency domain of nlbpsk.mat (top) and n2bpsk.mat (bottom) 
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Figure 7-10 - Mag(SCF) vs. alpha for nibpsk.mat (top) and n2bpsk.mat (bottom) 
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Figure 7-11 - Output of speccoa routine for nl bpsk.mat and n2bpsk.mat 


B. ARL:UT COLLECTED DATA 


The portion of the November 1995 Austin tests that are used in this work were 
collected by four stationary observers in the Austin area. The emitter was a push-to-talk 
narrow-band FM (NBFM), | W transmitter located near the center of the four observers. 
The GDOP of the scenario, given the positions of observation and the stationary target, 
was excellent at less then two. 

During the test, each observer records ~800 ms of data every 5 seconds. This data 
was stored to magnetic tape for post test analysis and further off-line testing. It is used to 
compare both TDOA determination and geolocation performance of the ARL:UT 
algorithms and those developed here. In order to best understand the output of ssca, 
plotscf and speccoa for the ARL:UT data, it is necessary to follow the signals’ processing 


from antenna to storage media. The processing at each observer is identical with respect 


7\ 


to archiving the data for later use and thus a single explanation serves for all four 


observers in this case. 
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Figure 7-12 - Beginning of ARL:UT signal processing 


Following reception, the 10 kHz bandwidth signal is heterodyned to an 
intermediate frequency (IF) of 21.4 MHz. At the IF, the 10 Msamples/second A/D 
converter aliases the signal at +1.4 MHz and -1.4 MHz with additional spectra added at 
intervals of the sampling rate. The spectrum at 1.4 MHz is downconverted to baseband 
as the product of an appropriate complex exponential multiplication. Its mirror at -1.5 
MHz is shifted to -2.8 MHz in the process. This series of steps is summarized in Figure 
7-12 with a spectrum relatively representative of that of the test data. 

At this point, the data is decimated in time by a factor or 1024. This step again aliases the 
signal at intervals of the new sampling rate, f, = 9765.625 Hz. The result is a set of four 
spectra in the interval - f, /2 to + f, /2. Figure 7-13 shows the results of this process in the 
time domain and frequency domain. Note the relationship between the highest negative 


frequency spectral line at approximately -440 Hz and the highest positive frequency 


a2 


spectral line at approximately 4440 is exactly fs/2 as they are related, aliased versions of 
the spectrum. The lowest negative frequency at approximately -4460 Hz and lowest 
positive frequency at approximately 420 Hz share that same relationship for the same 
reason. The frequencies are approximated above as they change slightly throughout the 
test and vary by observer as well. 

The corresponding SCF has features along the cyclic frequency axis that 
correspond to the various combinations of separation between the four spectral lines in 
the PSD of the signal. The greatest contribution to the SCF occurs at a separation of f, /2 
given the two constant relationships between the spectral lines sited above. At a = f; /2, 
all four lines contribute to the SCF thus the a = f, /2 feature contains the most energy. 
This is best illustrated the SCF and the magnitude vs. cycle frequency plots in Figure 7- 
14. Note the greatest peak in the SCF occurs at f; /2. 

Finally, using this greatest feature, speccoa finds the TDOA between observers. 
While other features may also lead to reliable answers, the f; /2 feature is a constant 
throughout the testing. All observers have this feature and it is always the strongest 
feature of the SCF. An example of the output from speccoa appears in Figure 7-15. It 
depicts the TDOA between the master node and slave #1. Specific results from the testing 
including some additional steps taken for more resolution in the TDOA computations 


appear in Chapter IX. 
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Figure 7-13 - Time (top) and frequency (bottom) plots of master signal, time 409400 
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Figure 7-14 - SCF (top) and Mag(SCF) vs. alpha for master signal, time 409400 
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Figure 7-15 - Output of speccoa for master and slave #3, time 409510 
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VII. GEOLOCATION SOFTWARE WORKBENCH BLOCKS 
A. SIMULINK® BLOCK CONSTRUCTION 


In addition to the testing of the combination of SPECCOA and Rusu’s closed 
form TDOA geolocation solution, the algorithms serve as the first blocks readied for 
construction for a developing geolocation software workbench. The details of their 
construction are outlined below in theory. 

Blocks that may be used in simulations run in the Simulink® environment 
conform to a general, highly flexible standard creates by The Mathworks, Inc. 
programmers. They may be implemented in one of three ways: graphically, through m- 
files or through another language such as ANSI C which produces a mex-file. No matter 
the method, the result is an S-function that operates according to the rules and protocols 
of the Simulink® environment.[24] 

Graphical creation of an S-function involves connecting the existing blocks that 
perform the needed function, collecting them as a single entity and masking that 
collection as the new block. This procedure forces Simulink® to create an S-function 
that describes the new system. This is by far the easiest method of creation for a speccoa 
block and a closed form TDOA geolocation block. The other two methods, m-files and 
mex-files, are appropriate in some instances such as certain state-space systems or other 
applications where the system can easily be written as a set of differential equations. 
Though more complex to create, S-functions created as m-files and mex-files simulate 
considerably faster. However, with the addition of the Simulink Accelerator®, 


graphically created S-functions enjoy some improvement in execution time as well.[24] 


B. SPECCOA BLOCK 


All the functions needed to use the speccoa code written for the testing described 
in both Chapters VII and IX, exist in Simulink® or one of the toolboxes. The Signal 
Processing Toolbox ® and DSP Blockset® contain the necessary buffers and functions 


needed by the speccoa code. A user defined Matlab® function block serves as the 


Te 


vehicle for the speccoa code to be used in the Simulink® environment. Connecting the 
buffers with the user defined function block constitutes the majority of the speccoa block. 
In general, a geolocation workbench block of speccoa requires two scalar 
sequences, the time samples of the two signals and one scalar parameter, the cycle 
frequency of interest as inputs. Time domain samples feed a buffer of length defined by 
the user in order to achieve the desired resolution. These time domain samples come to 
the speccoa block from an appropriately designed data conversion block. Regardless of 
its input, the conversion block must provide the speccoa block with single real or 
complex values for each simulation time period. In addition, the characteristic cycle 
frequency enters the block following the operation of a characterizer (defined elsewhere) 
to complete the inputs. The output of the speccoa block is a 2 x 1 vector, a time-stamp 
anda TDOA. Figure 8-1 shows a schematic of the speccoa block. Its role in a possible 


simulation scenario appears in Figure 8-3. 
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Figure 8-1 - Schematic of speccoa geolocation workbench block 
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c: GEOLOCATION BLOCK 


The closed form TDOA geolocation block is the more complex of the two 
algorithm blocks. Its inputs include four 3 x | observer position vectors with x, y and z 
parameters for each observer. In addition, it must have the output of three speccoa 
blocks. The geolocation block’s output is a 4 x 1 vector that includes the target’s x, y and 
z coordinates and a time corresponding to that geolocation. Figure 8-2 illustrates a basic 
diagram of the block. Its role in a possible simulation scenario appears in Figure 8-3 as 


well. 
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Figure 8-2 - Closed form TDOA geolocation solution block schematics 





D. SIMULATION SCENARIO 


One possible scenario that might be simulated using the two blocks described 
above is the November 1995 Austin test of the ARL:UT system. Using the signal data 


files and position files for each of the four observers, both the speccoa block and the 


19. 


closed TDOA form geolocation block are testable. The addition of a characterizer is 
necessary to feed the speccoa block. This function could be performed by a user defined 
block that calculates the maximum value of the SCF for each value of cycle frequency as 
plotscf plotted (see Chapter VII). The maximum value of this function could be the 
characteristic cycle frequency used in the speccoa block. In addition, the geolocation 
block could feed a plotting routine to display each geolocation. A diagram of this 


implementation appears in Figure 8-3. 
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Figure 8-3 - Simulation scenario 


Clearly the Simulink® environment is ideal for testing the many aspects of 
geolocation from data processing to parameter computation to geolocation estimation. 
The user defined MATLAB® function block provides for easy introduction of m-file 
routines into the S-function domain. It allows the interchange of various blocks in any 
scenario. This fact alone lends great flexibility and power to a geolocation software 
workbench centered around this environment. Other blocks such as pre-filters for the 


data or logic to evaluate the quality of the TDOA(s) or geolocation can also be 
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constructed. In addition, data from many different sources 1s handled simply by 
modifying the data conversion block to accept a different form of input. Its output 
remains the same. Using the Simulink® environment eliminates the need to define 
standards for module construction and provides an established simulation engine with 
many options. It is surely the most expeditious and intelligible implementation of a 


geolocation simulation system. 
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IX. RESULTS 


A. BACKGROUND 


Two comparisons are made between the cyclostationary TDOA processor 
measurements and those of the ARL:UT system. Both the TDOA(s) and the resulting 
geolocations serve as benchmarks of performance. In comparing the two systems’ 
results, theoretical TDOA(s) calculated from the target’s position and the positions of the 
four observers, using a value of 2.998 x 10° for c, the speed of propagation, are the 
baseline for comparison of the TDOA measurements. The target’s actual position 
throughout the test acts as the baseline for comparison of the geolocations. 

The TDOA(s) measured by both systems are biased by the radios used in the test. 
This constant bias was discovered in post test analysis and removed from the ARL:UT 
published results. However, the biases remain intact here to simplify the calculations 
and comparison. The geolocations suffer from the same biases and are also left 
uncorrected for the comparison. 

The comparison consists of 25 time periods of data from the November 1995 
Austin test. Because all four observers recorded data during these times, geolocations in 
three dimensions are possible. The process of arriving at three TDOA(s) and a 
geolocation involves plotting the SCF(s) for the four observers, choosing a cycle 
frequency, feeding that data to speccoa and finally running the closed form solution for 
the geolocation. The plots for the four SCF(s), the four magnitude of SCF vs cycle 
frequency and the three speccoa outputs appear in Appendix B for each of the 25 time 
periods. As noted in Chapter VII, the feature of the SCF at f,/2 provides speccoa the 


needed cycle frequency in every time instance. A summary of the results appears below. 


B. TDOA COMPARISON 


The comparison of TDOA(s) between the master node and slaves 1, 2 and 3 
appear below in Figures 9-1 through 9-6 below. Figures 9-1, 9-2 and 9-3 merely plot the 


TDOA(s) in the order in which they were determined. 
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Figure 9-1 - TDOA(s) for master node and slave #1, ARL:UT (top), speccoa (bottom) 
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Figure 9-2 - TDOA(s) for master node and slave #2, ARL:UT (top), speccoa (bottom) 
85 


t+tettareye tt ye ttt4e¢+t+euntts+es 


TDOA (seconds} 





times 


TDOA (seconds} 





0 5 10 |S 20 29 
times 


Figure 9-3 - TDOA(s) for master node and slave #3 ARL:UT (top), speccoa (bottom) 
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Figures 9-4, 9-5 and 9-6 show the deviation in seconds, from the calculated, 
theoretical TDOA, for each of the 25 time periods along the horizontal axis. The top 
plots in each figure again show results from the ARL:UT test; the lower plots show 


results of the cyclostationary TDOA processor. 


Cc. GEOLOCATION COMPARISON 


While the algorithm in both instances is identical, a comparison of the results 
serves to further illustrate the performance of the two different TDOA determination 
algorithms. The method of comparison is merely the geolocation performance. 
Unfortunately, even given the wide variety of TDOA(s) determined by the 
cyclostationary TDOA processor, the version of the closed form TDOA geolocation 
solution coded in MATLAB® produced nearly the same result for all 25 time instances. 
While this consistency is commendable, the various TDOA(s) should have produced 
fixes separated by several kilometers. This lack of definite variation in the geolocations 
reflects poorly on the quality of the algorithm’s outputs. 

The actual location of the transmitter is given by the coordinates 
Xt = -741941.474 m, y; = -5462120.413 m and Z% = 3198242.718 m. The MATLAB® 
solution produces xy = -740591 m, ym = -5443692 m and zy = 3187478 m. The only 
change to these rounded integer values over the course of the 25 time periods occurs 1n 
the decimal places. Clearly this is a very poor fix in addition to the lack of variation. 
The ARL:UT system successfully found the coordinates to within 100 m on all three axis 
for all 25 time periods. More often than not, the errors in the ARL:UT calculated values 
were on the order of 10 to 20 meters. 

Because the MATLAB® algorithm did not perform properly, a valuable 
comparison is impossible. The TDOA(s) discussed above are the only viable measure of 


comparison available. 
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Figure 9-4 - Deviation from theoretical for TDOA between master node and slave #1, 


ARL:UT (top), speccoa (bottom) 
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Figure 9-5 - Deviation from theoretical for TDOA between master node and slave #2, 


ARL:UT (top), speccoa (bottom) 
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Figure 9-6 - Deviation from theoretical for TDOA between master node and slave #3, 


ARL:UT (top), speccoa (bottom) 
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D. COMMENTS 


Clearly the cyclostationary TDOA processor does not perform as predicted by 
theory and as verified by simulation. Due to the decimation of the original time domain 
signals, several additional processing steps are necessary in order to achieve the needed 
resolution in the TDOA computation. These steps include resampling the signal at a 
much higher rate in order to measure a TDOA on the order of 10°. This step Causes 
problems for the speccoa code. 

The average values for the two TDOA determination methods appear in Tables 
9-1 to 9-3 below. Also in the table are the values for the average deviation of the 
TDOA(s) from the theoretical values. 

The code for the geolocation algorithm obviously has flaws as well. Its inability 

. to distinguish between the widely varying sets of TDOA(s) with widely varying 


geolocations indicates a fundamental programming flaw. 


Algorithm TDOAm) | Theoretical TDOA,y; 
Cyclostationary TDOA -1.308e-5 1.434e-6 1.165e-5 
ARL:UT System -1.805e-6 1.434e-6 3.66¢e-7 







Table 9-1 - TDOAy) summary 


Algorithm TDOAy? | Theoretical TDOA)> 


Cyclostationary TDOA 8 .44e-7 3.828e-6 2.984e-6 


ARL:UT System 2.697e-6 3.828e-6 1.131e-6 






Table 9-2 - TDOAy2 summary 
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Algorithm TDOAy3 | Theoretical TDOAy;3 
Cyclostationary TDOA -5.869e-6 -3.124e-6 2.745e-6 
ARL:UT System -3.47le-6 3.124e-6 3.47e-7 


Table 9-3 - TDOAm3 summary 
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x CONCLUSIONS 


A. OBSERVATIONS 


The cyclostationary TDOA processor in combination with the closed form 
geolocation solution in MATLAB® does not outperform the CAF based system. 
Because the results deviate significantly from the theory, several additional tests are 
necessary to isolate the cause of the anomaly. First, speccoa must run with fully 
decimated ARL:UT data. Because the decimation exceeds the resolution needed to 
determine the TDOA(s), properly the result should be zero in every instance. It is not. 
The algorithms returns random, obviously unreliable answers for the real world test data. 

The next step is the verification of zero delay with generated data. The 
generated BPSK signals with zero delay produce zero as a result with speccoa. Thus, 
speccoa indeed functions properly even at delays below the resolution of the sampling 
rate. However, it does not function properly with the ARL:UT data. 

Finally, since the SPECCOA algorithm reduces to a GCC at @ = 0 in the absence 
of Doppler shift, speccoa can be tested as a GCC to reproduce the ARL:UT results. This 
test produces random, unreliable results yet again with ARL:UT data. It does, however, 
produce the correct results with the generated BPSK signals once again. 

To summarize, plotscf and the cycle frequencies it produces are verified by the 
correct plots of the generated signals and the correct determination of the dominant 
cyclic features of the BPSK signals even in the presence of significant noise. Thus, the 
cycle frequencies used in speccoa for both the generated signals and the ARL:UT data 
are deemed reliable. Speccoa on the other hand, produces perfect results for the 
generated signals at all significant cycle frequencies and at © = 0. However, it does not 
produce consistent results for the dominant cyclic feature in the ARL:UT data; nor does 
it produce results at other cycle frequencies that are consistent with the fs/2 feature’ s 
results. 

These additional details lead to one most likely conclusion. The ARL:UT data 
does not reach the algorithms in the correct form from storage to processing. Most 


probably, the format is incorrect. In its reading into the MATLAB® environment, the 


Ze 


data is essentially corrupted producing such seemingly contrary results. The generated 
signals are generated in the MATLAB® environment, producing no such format 
problems. The ARL:UT data is originally stored in UNIX format, byte reversed for PC 
usage and then loaded into the MATLAB® environment for processing. It is this last 
step that is believed to cause the anomalies in the results. 

As a positive outcome, the testing did provide an idea for the computational 
complexity that a cyclostationary TDOA system would require. The process from 
loading the signals to finding a geolocation required approximately 2 minutes per 
geolocation on a Pentium Pro’™ based machine with 64 MB EDO RAM. This is not 
tactically feasible. While coding the system in C would increase the speed significantly, 
the computations of the system would remain the limiting factor in its implementation. 

The cyclostationary ‘TDOA processor holds promise for near future 
implementation. The problems encountered are tractable. All that remains 1s the 


computational power to run the system in a tactical environment. 


B. AREAS OF FURTHER STUDY 


Determining the proper method with which to load the ARL:UT data into the 
MATLAB® environment remains the highest priority for solving the major problems of 
the current algorithms. If this is the cause of the significant inconsistency encountered, 
its solution will surely improve the results to accuracies predicted in theory. 

The polynomial fit used by the resample function in MATLAB® needs to be 
examined closely to determine the extent of phase error it introduces to the resampled 
signal. This would allow for the development of a more accurate method of 
interpolating back to a sampling rate that provides sufficient resolution for the TDOA 
measurement. 

The closed form solution error must be found and corrected. It is likely that the 
error involves the reintroduction of c = 2.998 x 10° once the contributions of the 
TDOA(s) are mapped to the three coordinates. In his solution, Rusu sets c = | for 
simplicity. While this works for the derivation, it most likely must be given its proper 


value before the final solution can be reached. 
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Once the two algorithms were corrected, using real world signals from a more 
hostile signal environment for comparison would highlight the benefits of 
cyclostationary processing best. 

Finally, further development of Geolocation Software Workbench blocks would 


accelerate the introduction of a simplified testing platform for all aspects of geolocation. 
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APPENDIX A. MATLAB CODE FOR ALGORITHMS 


The programs that follow constitute the major sections of code used in this work. 
They all contain headers and commenting for explanation purposes. They are described 


briefly in Table A-1 below. 


Loads headers and signal data, converts 4 bytes real, 4 














bytes imaginary to complex vector notation. 








Strip-spectral correlation analyzer. Calculates the SCF 





of the input sequence along diagonal lines according to 


the SSCA algorithm. 





Transforms the diagonal lines of ssca.m into lines along 






f and a. Plots both the SCF and the magnitude of SCF 






vs cycle frequency. 


Resamples time domain data to user specified value (p). 






Calculates the TDOA of two input sequences in terms 


of samples using SPECCOA algorithm. 






Takes positions of observers and three TDOA(s) to 






produce geolocation in WGS84 ECEF coordinates. 


Finds both solutions. User must eliminate one. 


Table A-1 - Summary of code 
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fy 8 6 2 A A 2 A A 2 2 Ae 2 ee 2 EE A A AE 2 OC Fe 2 he 2 he 2 2 2 a ie 2 eC I OE Ok 
Of 6 Fe 2 2 2 2 AE 2 ie 2 Ae Ae 2 2 ie 28 A 2 AR 2 AE A 8 OR ee 26 2 he 2 2A oe 2 ee 2 2 Ae 2 2 ie 2 6 2 2 6 2 2 2 EO eK 2k 


% 

% Name: David A Streight 

% 

% Naval Postgraduate School - Monterey, California 
% 

% \loaddata.m 

% 

% Operating Environment 

J Windows 95 

% 

% Description 

% Loads 32-bit float ARL:UT data into four header and data vectors. 
J Master site is obsl, other sites are 2, 3 and 4. 
% 

% Date of last revision 

J 20 December 1996 

% 

% Inputs 

% Four ARL:UT stored files 

Jo 

% Outputs 

% None 

% 


fy AE RR Ae I RRS hea 2 oe ate te Ae oI a ok rh ara lete  e t 


oy 28 28 2 2 2 28 28 AR eR RR A 2k oR Go ok AG 2 2G AS AS G8 OR OR OR OR OR RR 2 28 9 OR OR 9 Ae AG OR OF ARR AS OR ahs 0 OK 2S OAS 26 AS OF HR AS Ae ols ols IS IS oe ote 


clear all; 


% open four observer files read only 

obs! = fopen('F:\data\nov_ut~2\mcc\409400.rev’,'r'); 
obs2 = fopen(F:\data\nov_ut~2\berg\409400.rev','r'); 
obs3 = fopen(F:\data\nov_ut~2\beecaves\409400.rev','r'); 
obs4 = fopen(‘F:\data\nov_ut~2\expo\409400.rev’,'r'); 


% \oad 3 header variables into header vectors 
[header1] = fread(obs1,3,'float'); 
[header2] = fread(obs?2,3,'float’); 
[header3] = fread(obs3,3,'float'); 
[header4] = fread(obs4,3, float’); 


% load data into 4 data matrices col#1 - real col#2 - imaginary 
[x1] = fread(obs 1,[8192,2],'float'); 
[x2] = fread(obs2,[8192,2],'float’); 
[x3] = fread(obs3,[8192,2], float’); 
[x4] = fread(obs4,[8192,2],'float'); 
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% convert data matrices to complex vectors 
S] = x1(:,1)4+)*x1¢,2); 
B= X2(:, 1 )4+)*x2(,2); 
Be = X3(:,1)4+)*x3C,2); 
S4 = x4(:, 1)+)*x4(:,2); 


% clear xi's 
Clear x |: 
clear x2; 
clear x3; 
clear x4; 


% close all files 
fclose(‘all'); 
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fo 78 28 Ae AE AR ee Oe OR RA 2 2 ie 2 2 8 IR 8 OR RE 8 26 2G ig 26 Eke 2c 2 2 2 2 ke kk 2k 2 2 2 2k 2 2 2 2 ok i 2 2 ok 2 ok ok OK 2 2k 2 2 


oy 7 AE AR AR AE AE AE ER OE 8 2 ER OR 28 RS 28 ES 2 gg 2k 2 2 I 2 2 2 2 OR Ok 2 i OE OO OOK 


%o 

% Name: David A Streight 

%o 

% Naval Postgraduate School - Monterey, California 
Jo 

Yo ssca.m 

%o 

% Operating Environment 

% Windows 95 

% 

% Description 

% Calculates the auto spectral density function using the strip 
% spectral correlation algorithm (ssca) 

%o 

% Date of last revision 

% 18 November 1996 

Jo 

% Inputs 

% x - SO! 

%o 

% Outputs 

% SDF - the results of the ssca computation. Must be converted 
Jo and plotted by plotsdf.m 

%o 


AE A A Rohe Ae Ae Ae Re Ae A Ae Ae eA Re eh ASR AER Ke AAs te Re Ae A Re Ria ae ae ea ate ata Re tae ee 


py 28 A 6 A OO AE OO OO ee EK A EO 2 OR 2 Ae AS OE IG EC 2 2 Oe 2G 2 2G 2 2 A AGE Ae AG Ee AG 2G AE 2G 2A oe 2S oO AE OE 
R= isl 272. 


% prepare variables 


S = s1ze(x); 

N= sCl): % number of points in vector 

Np 8: % course FFT number 

xp = [x zeros(1,Np)]; % zero pad input variable for Np length slices 
w = (hamming(Np))’; % hamming window for course FFT 

W = (hamming(N))’; % hamming window for P point FFT at the end 
L=Np/4; % decimation factor 

P = floor(N/L); 


p = zeros(1,floor(P/L)); 
xT = zeros(1,Np); 

k = -Np/2:Np/2-1; 

dm = zeros(1,Np); 

XTi = zeros(N,Np); 
SXT = zeros(Np,N); 
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% find complex demodulates 


for p = 0:P-1 
xT = fft((xp(1,p*L+1:p*L+Np) .* w), Np); % take Np point FFT of windowed data 
roe x 1528) x PCL, 1:4) ]; - % fold fft to account for neg freq's 
aim = exp(-|*2*pi*k*p*L/Np); % calculate freq shifts for this set 
XTi(p*L+1,:) = xT .* dm; % freq shift and populate odd col's of interp 
matrix 
end 


% interpolate ***works for Np = 8 only at this point*** 
for v = 1:P-1 

XTi(L*v,:) = (XTi(L*v-1,:) + XTi(L*v+1,:))./2;  % fill in even col's linearly 
end 


Pen) = X Ti1(1,:): % hold next to last col for last col values 
Zofinish up 

pam | X°X°XO KS XGXXSX] : % matrix of original signal 

JWm = [W3W;W;:W;W;W;W;W];: % matrix of column-wise hamming 
window 

SXT = fft(XTi .* xm, N); % take N-point FFT of products 

20 =SXT.; % transpose after FFT 

sl = (SDF. 1ePaliSDFC-,P:N)]; % fold fft to account for negative freq's 
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Oy FE A AE A A A 2 2 2 AS 2 2 CAR HR I AS OR 2h 2 AR AR OR AS 2 EE 2 AC 2 CE 3 2 2 2S 2 gO ig 2g 2 2 So aK ok OK 2 2 2k 


Oy 7 AE AE 2 eg 2 2 2k ee 2 AC OE 26 AS 2S 2 2S iS 2S 28 AE 2 2 AS AS DS 8 24g IC 2s 2s fe 2S 2 fe 2 AS 2S 2 fe 2g 2 2 ES 2 2 og 2 8 a A 2 a OE 2 2 OK OE 2k 


% 
% Name: David A Streight 


% 

% Naval Postgraduate School - Monterey, California 
% 

% PLOTSDF: Convert and plot SSCA output 

%o 


% Operating Environment 
% Windows 95 

% 

% Description 

%o 

% 

% 

% 

% Date of last revision 

- % 15 November 1996 

Jo 

% Inputs 

% Matrix of Spectral Density Function values from the SSCA algorithm 
% SDF 

Jo 

% All variables used in the SSCA computation (no clear executed) 

Jo 

% Outputs 

% Plot in 3D of the SDF function drawn along lines of cycle frequency 
% 


Of A 2 2 2 Ae ee ee a ae ie 2 oe 9 2 ae Ae 2 2 AS A hs Aa eR te 2 ta a Eats 2 Aa Aa os a ee 


Oy 7 EE AG AE OE 2 2 AG ie 2 Fa 2 OE ae A AG Oe AG Ae AE 2K 9 2 ie A OD 2 OF. DO a ae Ae A ae he AGI eke tao oRe ate he 2 oko ota 2h Pe ofa ee ok a 


% determine the span of alpha 
alphamax = 2*N - 2*N/Np - 2; 
alphamin = -2*N; 

alphaspan = alphamax - alphamin; 


% setup intermediate matrices 
1 = zeros(alphaspan+1,Np); 

X = zeros(alphaspan+1,Np); 
Y = zeros(alphaspan+1,Np); 


%o 
maxalpha = zeros(1],alphamax+1); 


% populate matrices from SDF 
for q = -N/2:N/2-1 
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for k = -Np/2:Np/2 - 1 
[(2*N*k/Np+2*q+2*N+1,k+Np/2+1) = abs(SDF(k+Np/2+1,q+N/2+1)); 
end 

end 


% normalize magnitudes 
{ = I/max(max(])); 


% populate f and alpha coordinate matrices 

= 11:8; 

y = ones(1,8); 

for 1= l:alphaspan+] 
X(i,:) = 2*N*(j-Np/2-1)/Np-(i-2*N-1)/2; 
nea(ies) —((1-2* N-1).*y)/2; 


end 


% plot the positive alpha half of the plane 

figure(1); 

surf(X(alphaspan-alphamax:alphaspan+1,:), Y(alphaspan- 
alphamax:alphaspan+1,:),[(alphaspan-alphamax:alphaspan-+1,:)); 
view(-105,60),grid,axis({-20000 10000 0 8000 0 1 )) 


% now plot alpha vs mag of SDF for those values 
for | = l:alphamax+1 

maxa(l) = max(I(alphaspan-alphamax+1-],:)); 
end 


maxalphas = nonzeros(maxa); 
maxalphas = maxalphas’; 

SZ = size(maxalphas); 
alphanum = 0:sz(1,2)-1; 
figure(2); 
plot(alphanum,maxalphas) 


clear 1; 
Clear J; 
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To 26 EAE EK A AR AER AR AR KOR 2 A A OK 8 2 2 2K 2k 2 2 2 2 2 2 RR 2 OE 2 2 2 2K 2K OK 2 2 2 2 OK ok KOK Ok 
Oo 7 7 AE AE AR AR AR EE AE EA Og A OR 2 A Og GE OS 2 OG OK 2 2 2 oe Oe 2 2 2k 2 ge 2 OR kK OK KO 2 OK KOK 


%o 

% Name: David A Streight 

Yo 

% Naval Postgraduate School - Monterey, California 

%o 

% SPECCOA - Cyclostationary Time Difference of Arrival Algorithm 
Yo 

% Operating Environment 

% Windows 95 

Jo 

% Description 

% Produces an estimate of the time difference of arrival of a 
% signal of interest between two spatially separated receivers 
% in terms of multiples of the sampling frequency 

% 

% Date of last revision 

% 06 January 1997 

%o 

% Inputs 

% p-up sample integer 

% q- down sample integer 

% terms - number of terms used in linear interpolation 

% N-number of samples of decimated sequence to use 

% a-cycle frequency of interest in terms of decimated fs/N 
% SOJI - signal captured at receiver | (N samples) 

% SOJ2 - signal captured at receiver 2 (N samples) 


%o *both SOI(s) sampled at rate fs samples per sec 
%o 
% Outputs 


% Plot of SPECCOA function versus samples. Peak is TDOA estimate 
%o 


Wo 78 78 AE AR Ae AR A AE AR 2 8 A a KK OK OK OK OK 2 2 KE OK ok KK OK OK og 2 ek og og og 8 gO OK OK OK OK OK OB 2 OBE 8 OK OK OK KKK 
Oo 78 FE A 28 AE RK 2 8 2K OK 2 OK OK RG A ge ie 2 2 2 2 2 2 8 ok OO og 8 8 2 2 2 ie ie gs 2g Ok 2 2g 2 OK OK Ok OK OK OK OK 


loaddata; % load decimated signals 


% resample variables 


p = 256; % up sample 

=I; % down sample 

terms = 3; % terms in the linear interpolation 
I = p/q; % total interpolation 


% cyclostationary variables 
a = 2049 * J; % cycle frequency at new sample rate 
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N = 4096 * I; Yo length of sequence at new sample rate 
fs = 9765.625; 


% prepare loop variables 


SXT = zeros(1,N); %o SCF of SOT] at alphaO 
SYXT = zeros(1,N); % Cross SCF at alphaO 
Sig] =s$1(1:(N/I),1).; % set up vectors 


Sig2 = s2(1:(N/D,1).’; 


SOI] = resample(Sigl,p,g,terms); | _% resample and set as SOI] 
SOI2 = resample(Sig2,p,q,terms); | _% resample and set as SOI2 


Xf = fft(SOI1,N); e % freq spectrum of SOI] 
XT = fftshift(Xf); % adjust for MATLAB fft 
Yf = fft(SOI2,N); % freq spectrum of SOI2 
nel = fftshitt(Yf); % adjust as above 


XT = [zeros(1,a/2) XT zeros(1,a/2)]; % zero pad at each end for correlation loop 
YT = [zeros(1,a/2) YT zeros(1,a/2)]; % ditto 


SXT(1,:) = XT(1,1:N) .* conj(XT(1,a+1:N+a)); 
SYXT(1,:) = YT(1,1:N) .* conj(XT(1,a+1:N+a)); 


isp = ifft((conj(SXT) .* SYXT), N); % IFFT of complex product from above 


s = -N/2:N/2-1; % time in samples 
Kemel= Expy pita*s/N); % freq kernel at alphaO 
sp = fftshift(isp) .* kernel; % adjust and multiply 
plot(s,real(sp)) % plot to find TDOA 
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Of 7 78 2 AR 28 AR A 2 ee RE 8 gE OR OR RE 8 Ek 8 2 8 he 2 2 ie AS eG kG 2S 2 2 Oe 8 8 OR OK 8 OK OKO 6 2K 


Op 7 28 2 AR 28 Ag A 2 ge ie he EE se AE OR 8 8 AE eG 8 8 8 ie A 2 gg ke A 8 2 OO 2 2 2k 6 OB 6 Oe OK OK oe OK 


Jo 

% Name: David A Streight 

Jo 

% Naval Postgraduate School - Monterey, California 

Jo 

% Closed Form TDOA Geolocation Solution 

Jo 

% Operating Environment 

% Windows 95 

Jo 

% Description 

% Produces two WGS84 based 3-dimensional geolocation solutions for a transmitter 
% given the location of each of four observers and the TDOAs for their observations 
% Data outside the algorithm must be used to eliminate the ambiguity. 
Yo 

% Date of last revision 

% 10 October 1996 

Jo 

% Inputs 

Jo WGS84 based x, y and z coordinates for each of four observers 

% xi, yi zi,1= 1,2,3,4 

Jo 

% TDOAs for the observers 

%  Dy,1=1,2,3,4 andj = 1,2,3,4 where 1! j and Diy =-Dyji 

Jo 

% Outputs 

Jo WGS84 based x, y and z coordinates for the transmitter 

%  rtl and rt2 (vectors) 

Jo 


A ois A a As he AS AG Ae AS Ae Ae AS ots “An oA oh eons Ae Ae 2A obs An Ae eke Aah Aa oe 2h Ae he Pe Ae Ae Ae Ra te ee 


» Gy FH FA AA ARR 8 OR AE OR AE AB AR AR RoR AS AS 8 OB OR 2h 2 AS OR aie AS AS 2s he AS AS 8 ae 2 AS AS hs AS AS AS oe Se oBe he ee che ote he che che obs hv cde Pain ote chalet 


% positions 

xt = -741941.47; yt = -5462120.41; zt = 3198242.728; 
X1 = -741203.52; yl = -5456323.70; z1 = 3208396.55; 
X2 = -735740.31; y2 = -5467456.78; z2 = 3190514.57; 
X3 = -754177.92; y3 = -5459066.31; z3 = 3200764.99; 
X4 = -731248.92; y4 = -5463237.87; z4 = 3198800.64; 


% set up inputs and utility vectors 


tl =0; % relative times of arrival using observer | as 
baseline 

{2 =—-Di2° % D12 is defined as tl] - t2 thus t2 =-D12 if tl =O 
t3 = -D13; % ditto 
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t4=-D14; % ditto 


rl =([x1;y1;z1]; % observer position vectors 
ea [X2-y 2-72 |: 
fe (X55 yoeZ a: 
r4 = [x4;y4;z4]; 


Euchid_sq] =x142 + y142+2z1%2; % squares of Euclidean distances 
Epoicesg. = x2°2 4+ y2"2 + 22% 2: 
ucidisg3s = x3°2 + y3"2 + z3“2; 
Euclid_sq4 = x4%2 + y4%2 + 242: 


% set up solution matrices 
A = [x1-x2 yl-y2 z1-22;x2-x3 y2-y3 22-z3;x3-x4 y3-y4 z3-z4]; 


q = [t1-t2;t2-t3;t3-t4]; 


Se= zeros(3,1); 

Slee = (Buchdesq! - Euclid sq2 - tl*2 + t2“2).*0.5; 
S(ze)) —(uchid sq2 - Euclid sq3 - t2“%2 + t3*2).*0.5; 
s(3,1) = (Euclid_sq3 - Euclid_sq4 - t342 + t4%2).*0.5; 


g = A\q; % use matrix division (MATLAB)... 
% in lieu of the inverse function... 
h = A\s; % to solve this part 


% find terms for at42 + bt +c 
a = ((g(1,1))42 + (¢(2,1))42 + (2(3,1))*2) - 1;% find quadratic term 


pee, * hh) +(e. * rl) + tl; % find first order term 
c= (norm(h = rl).*2) - t142; % find constant 
0) = |B evel p % form polynomial in t 


% find the two times as roots of this equation 
t = roots(p); % find roots of at42 + 2bt +c 


% substitute into t-parametrized equation for r and solve for real part of two solutions 


nole—wreate tc leh) + h); 
Ft2 = teal(e .* ((2,1|) + h); 
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APPENDIX B. PLOTS FOR ARL:UT TEST DATA 


The following 150 figures are the result of testing 25 the cyclostationary TDOA 
processor for 25 time periods. For each time period, 6 figures appear: one SCF plot for 
each of the four observers grouped two to a figure, one magnitude of SCF vs cycle 
frequency for each of the four observers grouped two to a figure and three SPECCOA 
plots grouped two and one to a figure. 

For the SCF plots, only the positive cycle frequency values are plotted. The 
spectral frequency axis runs from the upper left to the lower right; the cycle frequency 
axis runs from the upper right to the lower left. Both the spectral frequency and the 
cycle frequency are plotted in terms of f,/N where f; = 9765.625 Hz and N = 4096. The 
magnitude is of course the axis pointing to the top of the page. The magnitude of the 
plots is nomalized by the maximum value of the SCF and thus ranges from 0 to I. 

The magnitude of the SCF versus cycle frequency plots also display only the 
positive cycle frequencies. The horizontal axis is cycle frequency in terms of f,/N as 
above. The vertical axis is again normalized magnitude. These plots merely represent 
looking at the SCF plots from the cycle frequency axis side. They aided in choosing the 
best cycle frequency for SPECCOA. 

Finally, the SPECCOA plots show the output of the speccoa code. The 
horizontal axis represents the TDOA in terms of the numbers of samples of delay. The 


horizontal axis is once again magnitude. 
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Figure B-1 - SCF for master (top) and slave #1 (bottom), time 409400 
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Figure B-2 - SCF for slave #2 (top) and slave #3 (bottom), time 409400 
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Figure B-3 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409400 
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Figure B-4 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409400 
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Figure B-5 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409400 
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Figure B-6 - SPECCOA for master and slave #3, time 409400 
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Figure B-7 - SCF for master (top) and slave #1 (bottom), time 409405 
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Figure B-8 - SCF for slave #2 (top) and slave #3 (bottom), time 409405 
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Figure B-9 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409405 
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Figure B-10 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409405 
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Figure B-1 1 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409405 
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Figure B-12 - SPECCOA for master and slave #3, time 409405 
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Figure B-13 - SCF for master (top) and slave #1 (bottom), time 409410 
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Figure B-14 - SCF for slave #2 (top) and slave #3 (bottom), time 409410 
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Figure B-15 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409410 
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Figure B-16 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409410 
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Figure B-17 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409410 
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Figure B-18 - SPECCOA for master and slave #3, time 409410 
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Figure B-19 
-19 - SCF for 
master (top) and slave #1 (bottom) 
m), time 409415 
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Figure B-20 
-20 - SCF for sl 
ave #2 (top) and slave #3 (bottom) 
m), time 409415 
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Figure B-21 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409415 
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Figure B-22 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409415 
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Figure B-23 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409415 
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Figure B-24 - SPECCOA for master and slave #3, time 409415 
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Figure B-25 - SCF for master (top) and slave #1 (bottom), time 409420 
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Figure B-26 - SCF for slave #2 (top) and slave #3 (bottom), time 409420 
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Figure B-27 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409420 
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Figure B-28 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409420 
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Figure B-29 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409420 
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Figure B-30 - SPECCOA for master and slave #3, time 409420 
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Figure B-31 - SCF for master (top) and slave #1 (bottom), time 409425 
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Figure B-32 
-32 - SCF for sl 
ave #2 (top) and slave #3 (bottom) 
om), time 409425 
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Figure B-33 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409425 


146 





0 900 1000 1500 2000 2500 3000 3500 4000 





0 500 1000 1500 2000 2500 3000 3500 4000 


Figure B-34 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409425 
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Figure B-35 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409425 
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Figure B-36 - SPECCOA for master and slave #3, time 409425 
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Figure B-37 - SCF for master (top) and slave #1 (bottom), time 409430 
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Figure B-38 - SCF for slave #2 (top) and slave #3 (bottom), time 409430 
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Figure B-39 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409430 
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Figure B-40 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409430 
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Figure B-41 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409430 
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Figure B-42 - SPECCOA for master and slave #3, time 409430 
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Figure B-4 
-43 - SCF for 
fac 
aster (top) and slave #1 (bott 
om), time 4094 
35 
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Figure B-44 - SCF for slave #2 (top) and slave #3 (bottom), time 409435 


pe 


Oo 


0.8 


O./ 


0.6 


Wo 


04 


Ors 


0.2 





0.15 


| | : 


} 


0 900 1000 1500 2000 2500 3000 3500 4000 


0.9 
0.8 
0.7 
0.6 
0.5 
04 
Os 
0.2 





Jf 
wren ee 
0 900 1000 1500 2000 2500 3000 3500 4000 


Figure B-45 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409435 
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Figure B-46 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409435 
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Figure B-47 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409435 
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Figure B-48 - SPECCOA for master and slave #3, time 409435 
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Figure B-49 - SCF for master (top) and slave #1 (bottom), time 409455 
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Figure B-50 - SCF for slave #2 (top) and slave #3 (bottom), time 409455 
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Figure B-51 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409455 
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Figure B-52 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409455 
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Figure B-53 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409455 
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Figure B-54 - SPECCOA for master and slave #3, time 409455 
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Figure B-55 
-55 - SCF fo 
r master (top) and slave #1 (bottom), t1 
, time 409460 
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Figure B-56 - SCF for slave #2 (top) and slave #3 (bottom), time 409460 
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Figure B-57 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409460 
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Figure B-58 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409460 
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Figure B-59 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409460 
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Figure B-60 - SPECCOA for master and slave #3, time 409460 
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Figure B-61 
-61 - SCF for 
master (top) and slave #1 (bottom), ti 
, ume 409465 
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Figure B-62 - SCF for slave #2 (top) and slave #3 (bottom), time 409465 
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Figure B-63 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409465 
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Figure B-64 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409465 


177 





-2 
-8 -6 4 -2 0 Z 4 6 8 


4 
x 10 
Figure B-65 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409465 
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Figure B-66 - SPECCOA for master and slave #3, time 409465 
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Figure B-67 - SCF for master (top) and slave #1 (bottom), time 409470 
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Figure B-68 - SCF for slave #2 (top) and slave #3 (bottom), time 409470 
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Figure B-69 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409470 
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Figure B-70 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409470 
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Figure B-71 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409470 
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Figure B-72 - SPECCOA for master and slave #3, time 409470 
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Figure B-73 - SCF for master (top) and slave #1 (bottom), time 409475 
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Figure B-74 - SCF for slave #2 (top) and slave #3 (bottom), time 409475 
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Figure B-75 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409475 
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Figure B-76 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409475 
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Figure B-77 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409475 
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Figure B-78 - SPECCOA for master and slave #3, time 409475 
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Figure B-79 - SCF for master (top) and slave #1 (bottom), time 409480 
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Figure B-80 - SCF for slave #2 (top) and slave #3 (bottom), time 409480 
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Figure B-81 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409480 
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Figure B-82 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409480 
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Figure B-83 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409480 
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Figure B-84 - SPECCOA for master and slave #3, time 409480 
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Figure B-85 - SCF for master (top) and slave #1 (bottom), time 409485 
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Figure B-86 - SCF for slave #2 (top) and slave #3 (bottom), time 409485 
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Figure B-87 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409485 
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Figure B-88 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409485 
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Figure B-89 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409485 
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Figure B-90 - SPECCOA for master and slave #3, time 409485 
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Figure B-91 - SCF for master (top) and slave #1 (bottom), time 409490 
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Figure B-92 - SCF for slave #2 (top) and slave #3 (bottom), time 409490 
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Figure B-93 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409490 
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Figure B-94 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409490 
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Figure B-95 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409490 
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Figure B-96 - SPECCOA for master and slave #3, time 409490 
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Figure B-97 - SCF for master (top) and slave #1 (bottom), time 409495 
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Figure B-98 - SCF for slave #2 (top) and slave #3 (bottom), time 409495 
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Figure B-99 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409495 
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Figure B-100 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409495 
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Figure B-101 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409495 
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Figure B-102 - SPECCOA for master and slave #3, time 409495 
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Figure B-103 
- - SCF for 
master (top) and slave #1 (bottom) 
m), time 409500 
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Figure B-104 - SCF for slave #2 (top) and slave #3 (bottom), time 409500 
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Figure B-105 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409500 
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Figure B-106 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409500 
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Figure B-107 - SPECCOA for master and slave #1! (top) and master and slave #2 (bottom), time 409505 





=. -6 4 -2 0 Z 4 6 8 


Figure B-108 - SPECCOA for master and slave #3, time 409505 
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Figure B-109 - SCF for master (top) and slave #1 (bottom), time 409505 
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Figure B-110 - SCF for slave #2 (top) and slave #3 (bottom), time 409505 
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Figure B-111 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409505 


224 





0 900 1000 1500 2000 2500 3000 3500 4000 





0 500 1000 1500 2000 2500 3000 3500 4000 


Figure B-112 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409505 
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Figure B-113 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409505 
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Figure B-114 - SPECCOA for master and slave #3, time 409505 
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Figure B-115 - SCF for master (top) and slave #1 (bottom), time 409510 
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Figure B-116 - SCF for slave #2 (top) and slave #3 (bottom), time 409510 
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Figure B-117 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409510 
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Figure B-118 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409510 
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Figure B-119 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409510 
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Figure B-120 - SPECCOA for master and slave #3, time 409510 
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Figure B-121 - SCF for master (top) and slave #1 (bottom), time 409515 
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Figure B-122 - SCF for slave #2 (top) and slave #3 (bottom), time 409515 
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Figure B-123 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409515 
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Figure B-124 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409515 
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Figure B-125 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409515 
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Figure B-126 - SPECCOA for master and slave #3, time 409515 
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Figure B-127 - SCF for master (top) and slave #1 (bottom), time 409520 
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Figure B-128 - SCF for slave #2 (top) and slave #3 (bottom), time 409520 
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Figure B-129 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409520 
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Figure B-130 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409520 
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Figure B-131 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409520 
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Figure B-132 - SPECCOA for master and slave #3, time 409520 
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Figure B-133 - SCF for master (top) and slave #1 (bottom), time 409525 
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Figure B-134 - SCF for slave #2 (top) and slave #3 (bottom), time 409525 
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Figure B-135 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409525 
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Figure B-136 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409525 
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Figure B-137 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409525 
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Figure B-138 - SPECCOA for master and slave #3, time 409525 
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Figure B-139 - SCF for master (top) and slave #1 (bottom), time 409530 
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Figure B-140 - SCF for slave #2 (top) and slave #3 (bottom), ttme 409530 
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Figure B-141 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409530 
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Figure B-142 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409530 
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Figure B-143 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409530 
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Figure B-144 - SPECCOA for master and slave #3, time 409530 
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Figure B-145 - SCF for master (top) and slave #1 (bottom), time 409535 
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Figure B-146 - SCF for slave #2 (top) and slave #3 (bottom), time 409535 
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Figure B-147 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409535 
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Figure B-148 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409535 
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Figure B-149 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409535 
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Figure B-150 - SPECCOA for master and slave #3, time 409535 
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